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Synopsis 


Kinetic simulations based on the Boltzmann equation with Fokker-Planck collision 
term have been quite useful in modeling multispecies plasmas. Numerical simula- 
tions of fusion plasmas often deal with many ion species at different temperatures 
and also the electrons whose characteristic velocity and time scales are very dif- 
ferent. Further, highly peaked non-Maxwellian distributions which arise out of 
auxiliary heating and reaction products are also encountered in the simulation of 
fusion plasmas. Such physical scenarios naturally add to complications in the nu- 
merical modeling process. Therefore it is desirable to represent different species 
including electrons on their cheiracetristic velocity scales to facilitate the study of 
energy transfer, auxiliary heating, interaction with the fusion reaction products and 
of other physical situations of practical interest. 

Fokker-Planck equations are solved to obtain the distribution functions, in gen- 
eral, by the finite difference method. Relatively less work has been done using finite 
element method for solving Fokker-Planck equations. Further, modeling species 
with highly peaked initial distributions, their subsequent interaction with other 
species and their evolution can be effectively dealt with by adapting the numerical 
grid. This leads to improvements in the accuracy of solutions, and in the calculation 
of such pareuneters as reactivities, etc. 

Some of the contributions of the present work can be listed £is given below: 

1. Multi-velocity and time scale Fokker-Planck equations in 2-D velocity space 
{v, $), where v is the speed and 6 the pitch angle, are introduced. These equa- 
tions are derived by using individual velocity cut-offs for ions and electrons. 
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Therefore electrons can also be treated together with ions as a general species 
with less number of grid points. 

2. Multispecies Fokker-Planck equations both in 1-D and 2-D are solved by 
Galerkin finite element method using conventional finite elements. Relative 
performance of various finite elements is studied. The Fokker-Planck equa- 
tions are also solved by the finite difference method. 

3. A method is presented to redistribute the nodes optimally on the numerical 
grid. Redistribution is based on the local discretization errors. 

4. The Fokker-Planck collision operators which have singularities at the bound- 
aries are approximated analytically. 

5. The Fokker-Planck models are used to calculate reactivities < crv > accurately 
during fast ion heating. 

detailed outline of the present work is given in the following paragraphs. 

Chapter 1 presents a review of literature on the numerical solution of Fokker- 
anck equations. A significant amount of work reported in literature is on the 
ite difference method rather than on the finite element approach. A brief review 
the area of optimal node distribution is also presented. Further, some studies 
rtaining to the < <tv > calculations and source and loss mechanisms that have 
en considered in the simulations are also reviewed briefly. This chapter ends with 
. outline of the work presented in subsequent chapters. 

Chapter 2 contains the governing equations and approximation of collision oper- 
ors at the boundaries. Multispecies Fokker-Planck equations in (u, 9) co-ordinates 
e written in their general form. These equations are recast into a dimensionless 
rm by using velocity and time constants which are chosen based on the charac- 
ristic velocity and time scales of different species. This formulation facilitates 
e treatment of electrons as a general species with the same computational effort, 
case of MaxwelHan background species, the Fokker-Planck coefficients are pre- 
nted in their simplified forms. The singularities present in the collision operators 
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at the boundaries (u = 0, 0 = 0 and tt ) are removed analytically. One-dimensional 
Fokker-Planck equations as derived from the 2-D equations are also presented. 

Chapter 3 is concerned with numerical solutions of multispecies 1-D Fokker- 
Planck equations. The equations are first solved by the Galerkin finite element 
method with the standard isoparametric finite elements namely, the linear, quadratic 
and cubic elements. Time stepping is done by the one step 0-methods (0 = 1/2 
- Crank-Nicolson method, 0 = 1- Backward implicit method). Computational 
aspects pertaining to the Fokker-Planck coefficients, use of higher order quadra- 
ture, velocity scale refinement, time step control strategy etc. are discussed. The 
spatial and temporal convergence of various finite elements and the time stepping 
schemes are studied. The results show that the quadratic and cubic elements are 
better than the linear elements in terms of accuracy and convergence. Multispecies 
calculations include treating electrons as a general species. The advantage of using 
multi- velocity scales is emphasized. Further, the 1-D Fokker-Planck equations are 
solved by the finite difference method. Their solutions are compared with finite 
element solutions. The collision operator without singularity at u = 0 is also im- 
plemented. The results indicate that the analytical approximation improves the 
accuracy of the solutions at the boundary u = 0. 

Chapter 4 deals with 1-D finite difference approach mainly from the view point 
of optimal distribution of grid points. The node redistribution is intended to reduce 
the TnaYimnTn of local discretization error. Higher order approximations for the 
derivatives are also considered and compared with standard central finite difference 
discretization. The discretization errors are compared with the exact values for 
some initial distributions. It is found that the maximum errors for Maxwellian 
and Gaussian distributions can be reduced considerably by (in most cases by more 
than an order of magnitude) the redistribution of nodes. The grid redistribution is 
illustrated with the relaxation of a Gaussian distribution to Maxwellian at various 
time intervals. 

Chapter 5 is mainly on the numerical solution of Fokker-Planck equations in 2-D, 
(v, 0) space. First, the equations are solved by the Galerkin finite element method 
using five different types of quadrilateral elements, viz. linear, complete quadratic 
and cubic, and incomplete quadratic and cubic elements. Relative merits of these 
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!ments axe investigated. As in 1-D, the time stepping is done by the one step 
■methods. Secondly, for purpose of comparison and to study the implementation 
the analytical approximation at the boundaries, the finite difference approach 
used. For this case Alternating Direction Implicit time stepping scheme is also 
eluded. 

Chapter 6 consists of some studies related to the simulation of D-T fusion plas- 
as. Emphasis is on the calculation of < cru > accurately under non-Maxwellian 
nditions arising out of fast ion heating scenarios. An improved reaction loss model 
rich correctly calculates the particle losses due to fusion reactions at different 
ergies is implemented. Loss and source mechanisms considered in the simulation 
dude reacting particle losses and the reaction source for a’s. The transfer of en- 
5y from a particles to background plasma is studied taking into account the a-a 
;eraction. 

Finally, Chapter 7 summarizes the results presented in the previous chapters 
d ends with some observations and suggestions for the extension of the present 
»rk. 
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Chapter 1 


Introduction 


1.1 Motivation and Objectives 


Fusion Plasmas are generally not in a state of local thermodynamic equilibritun. 
In such cases the charged particle distribution functions are non-Maxwellian. For 
many applications, the actual knowledge of the distribution functions is required. 
Examples include, among others, the study of auxiliary heating by neutral beams 
or microwaves, thermalization of alpha particles in D-T plasmas, and the study of 
runaway electrons in tokamaks. In all these cases, kinetic models based on the Boltz- 
mann equation with Fokker-Planck collision terms have been successful (Killeen, 
Kerbel, et ai, 1986). However, numerical simulations based on the kinetic equa- 
tions involve considerable computational effort due to the nonlinear nature of the 
governing partial integrodifferential equations. Further, widely different character- 
istic velocities and relaxation times encountered in plasmas together with highly 
peaked distributions arising out of auxiliary heating and reaction products add to 
the complications in the numerical modeling process. 

The finite difference method has been generally used for the numerical solu- 
tion of Fokker-Planck equations, with relatively less emphasis on the finite element 
method. In the present work both finite element and finite difference methods of 
discretization have been employed, with empheisis on the following aspects of the 
simulation process: 
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govornin^ FokkBr^Plsnck ©Qiisitions o-iid tliG RosGublutli pofcoiiti<ils ^rc 
cast in a dimensionless form using a cut-off velocity characteristic to each 
species. This allows the nodes for each species to be chosen independ(uitly of 
other species. 

ii. In putting the equations for a particular species in the dimonsionUvss form, 
the relaxation time of the species is used for normalizing the time variable. 
This facilitates each species to be integrated on its own characteristic time 
scale. This is possible in Fokker-Planck equations as the coupling between the 
different species is only through Rosenbluth potentials. 

iii. In spherical polar geometry, the Fokker-Planck operator has singularities at 
u = 0, and B — 0 and tt. These are eliminated analytically in the present 
work, and the resulting operator is discretized at the actual boundaries rather 
than at near-by points. 

iv. Multispecies Fokker-Planck equations both in 1-D and 2-D are solved by 
Galerkin finite element method. The relative performance of various finite 
elements is studied. Sufficient emphasis is put on the study of deasity and 
energy conservation, velocity scale refinement and its effect on the numerical 
oscillations, and tune step control based on the energy norm of the error in 
the computed solutions. 

V. Discretization errors over the velocity space are studied and a strategy to 
redistribute the nodes optimally to minimize the errors is presented. In view 
of the integrodifferential nature of the governing equations, the possibility 
of improving the accuracy using higher order interpolation and quadrature 
schemes is also studied. 

VI. A comparison of the Spitzer equipartition model with Fokker-Planck results, 
an improved fusion reaction loss term, and a few other aspects of energy 
transfer involving fast ions are also included. 

A brief outline of the work is presented in Section 1.4. In the next two sections 
the relevant literature for the present work is reviewed. 
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1.2 Numerical solution of Fokker-Planck equa- 
tions 

The literattire in the area of numerical solution of Fokker-Planck equations is re- 
viewed here. A stenting point could be the two books written 8 years ago (Killeen, 
Kerbel, et al, 1986; Dnestrovskii and Kostomarov, 1986) which give a wealth of 
literature and the history of the developments in the numerical simulation of fusion 
plasmas. These books provide a detailed description of several state-of-art Fokker- 
Planck computational models based on the finite difference method, together with 
many useful applications of these models to magnetically confined plasmas. 

The finite difference method of solution to Fokker-Planck equation is fairly estab- 
lished. Research has progressed step by step starting with single species calculations 
under simplifying assumptions to multidimensional models for multispecies. The 
Fokker-Planck equation was derived in spherical polar coordinates by Rosenbluth, 
MacDonald and Judd (1957) under the assumption that the distribution functions 
are azimuthally symmetric in (v, n, ip) space (/i = cos9, where 9 is the pitch angle). 
The potentials and the distribution functions are expanded in Legendre polynomi- 
als, which eventually lead to a system of coupled 1-D nonlinear integro-differential 
equations. Bing and Roberts (1961) solved the Fokker-Planck equation numerically 
for one ion species without taking into account the effect of electrons and other 
species. Their work also includes a feasibility study of using approximately sepa- 
rable solutions in (v, p) space. A similar work was attempted by BenDaniel and 
Allis (1962), with emphasis on approximating the distribution functions in a sep- 
arable form. Their study also involved modeling fusion applications for different 
parameters in mirror machines. 

The Fokker-Planck calculations in 1-D were later extended to more than one 
species. A study of energy transfer from hot ions to cold electrons was made by solv- 
ing Fokker-Planck equations for electrons keeping the ion distribution fixed in time 
(Killeen, Heckrotte and Boer, 1962). Steady state 1-D Fokker-Planck equations were 
solved both for ions and electrons by Fowler and Rankin (1966). The distribution 
functions are taken to be functions of energy. A time dependent Fokker-Planck code 
was later developed by Killeen and Futch (1968) to solve the 1-D kinetic equations 
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for ions and electrons in the velocity space v. The Fokker-Planck coefficients are 
extrapolated using the values at the previous time steps. The singularity at t; = 0 is 
avoided by solving the equations for the interior points only in the numerical grid. 
This code was used to model certain applications in mirror machines. 

Application oriented 1-D models were further developed by Killeen and Marx 
(1970). This work was also extended to 2-D in velocity space v and 6 for one ion 
species. The ion equation was solved by alternating direction implicit method. The 
coefficients in the 2-D Fokker-Planck equation were obtained by calculating the po- 
tentials g and h in their elliptic integral form. They have also discussed alternative 
ways of computing the potentials. Whitney (1970) solved the 2-D Fokker-Planck 
equation by expanding the solution in Legendre polynomials, which results in a sys- 
tem of 1-D equations. The equations are exact to all orders of Legendre polynomi- 
als. Until now the Fokker-Planck collision operator was expanded (nonconservative 
form) and then solved numerically. A new finite difference scheme which addresses 
mainly to the problem of distribution functions becoming negative was presented 
by Chang and Cooper (1970). They showed that their recommended scheme, based 
on conservative form of the collision operator, improves the accuracy of the solution 
and also possibly reduces the computational time. All the models/codes developed 

later on are mostly based on the conservative form of the Fokker-Planck collision 
operator. 

Full velocity space multispecies codes (Mirin, 1975; Killeen, Mrin and Rensink, 
1976) were developed to model appUcations in fusion plasmas for mirrors and toka- 
maks. These codes include many source and loss terms. The concept of hybrid 
models was also introduced in which the electrons are treated as isotropic. The 
2-D ion kinetic equations are solved either by the semi-implicit or by the alternat- 
ing direction implicit method, whereas 1-D equation is solved for electrons. These 
codes were also used to study the effect of anisotropic part of the ion-ion Rosen- 
bluth potentials on certain parameters of interest. Fowler, Smith and Rome (1978) 
developed a code (FIFPC) to study the slowing down process of fast ions in toka- 
m^s. This model treats the background plasma species to be Maxwellian. The 
Fokker-Planck equation is solved by the Crank-Nicolson and the alternating direc- 
tion implicit methods. A 2-D multispecies nonlinear Fokker-Planck package (FP- 
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PAC) was developed by McCoy, Mirin and Killeen (1981). In. this code derivatives 
in the velocity space are approximated by central differences. The colUsion operator 
which has singularities at the boundaries is also approximated numerically. Time 
stepping is made by semi-implicit or the alternating direction implicit or fully im- 
plicit methods. FPPAC is considered to be a standard 2-D Fokker-Planck code for 
plasmas in uniform magnetic fields. Later, they incorporated an improved finite 
difference scheme, which is similar to that suggested by Chang and Cooper (1970). 
According to them, FPPAC88, the newer version of the old code FPPAC, virtually 
eliminates the problem of distribution functions becoming negative (Mirin, Mccoy, 
et ai, 1988). 

Andrade and Oliphant (1981) solved 2-D Fokker-Planck equation by using a five 
point algorithm. They have calculated the coefficients with isotropic backgroimd 
species for all times. Spatial part is also included in their model. The resulting 
algebraic system of equations was solved by matrix factorization. For the 1-D 
time dependent Fokker-Planck equations several discretization schemes have been 
suggested by Larsen, Levermore, et al. (1985). According to them, their schemes 
possess certain advantages over the scheme of Chang and Cooper (1970). However 
their approach requires the coefficients of the equations in analytic form which can 
not be easily satisfied even in the case of 1-D multispecies models. 

Karney (1986) presented the numerical treatment of Fokker-Planck equation 
with particular reference to the study of rf driven current. He has discussed vari- 
ous time stepping schemes including the Chebyshev acceleration which is used to 
vary the time steps. Similar calculations for rf heating were made by O’Brien, Cox 
and Start (1986). It is suggested that unequally spaced mesh points in the finite 
difference or the finite element method may reduce the number of grid points re- 
quired. They ignored the singularities present in the collision operator by choosing 
the lower cut-off slightly away from the singular points or the boundaries. Recently 
McKenzie, O’Brien and Cox (1991) developed a code BANDIT3D based on the 
finite difference method to solve the 3-D Fokker-Planck equation for electrons in 
a tokamak. This work includes a Fokker-Planck model in (v, 0) together with one 
spatial variable v. Ions are assumed to be MaxwelUan. The equation is advanced 
in time using an operator splitting techmque in two stages. The code can be used 
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to study both kiuetic and transport physics. 

The progress in the area of finite element modeling of fusion plasmas has been 
slow till now. Fyfe, Weiser, et al. (1981) applied Rayleigh-Ritz finite element 
method to solve a linearized steady state Fokker-Planck equation in variational form. 
They used B-splines as basis functions for the phase space variables (energy and /j.), 
and solved the equation in a rectangular domain. A relativistic multiregion bounce 
averaged Fokker-Planck code for mirror plasmas was developed later (Matsuda and 
Stewart, 1986). This work includes solution of Fokker-Planck equation by Galerkin 
finite element method with specific applications to mirror plasmas. This code also 
uses one dimensional splines in energy and magnetic moment phase space. Poisson 
equations are solved directly by finite element method for the Rosenbluth potentials. 
Finite element method has also been used to model wave particle interactions for 


rf heating based on Fokker-Planck equations (Succi, Appert, et ai, 1986). 

Attempts have also been made to use the standard finite element library solvers 
to model plasma applications, mostly to solve a single kinetic equation. Shoucri, 
Fuchs and Bers (1987), Shoucri, Shkarofsky, et al. (1989) used standard 2-D finite 
element based partial differential equation solvers to model rf heating. These models 
mostly involve a fixed number of MaxweUians and seek the steady state solutions 


of Fokker-Planck equations for electrons. Harrison (1988) has applied moving finite 
element method to solve 1-D kinetic equations. Recently Gardner, Gardner and 
Zaki (1993) have modeled a few applications in which 1-D Fokker-Planck equation 
is coupled with other equations. The equations are solved by the Galerkin finite 
element method using the linear finite elements. 

Calculations of fusion reactivities < nt, > can be improved if the distribution 
funct. 0 ^ of the interacting species are known more accurately. The enhancement 
m the fusion reactivities due to auidUary heating was pointed out by Dawson 
^th and Tenney (1971). Since then fusion reactivity calculations have been per- 
orm^ under vanous conditions involving non-Maxwellian distributions (for exam- 
pe, utch Jr., Holdren, et at., 1972; Jassby and Towner, 1976; Slaughter, 1983). 
Over the ye^s, convenient, expressions and formulae have also been developed for 

Hain^^ cross-sections (Kvely, 1977; Bosch and 

ale, 1992). Faster methods have been developed to compute < nv >, which is a 
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5-D integral in velocity space, using Legendre polynomials (Cordey, Marx, et ai, 
1978). The 5-D integral can be simplified to 3-D integral if one of the distributions 
is isotropic (Meirx, Mirin, et al, 1976). 

Various studies have been made to investigate the effect of distortion in the ion 
distribution functions on the fusion reactivities. The improvement in the fusion 
reactivity and multiplication factor in the presence of fast ions has been studied 
recently by Niikura and Nagami (1988), Niikura, Nagami and Horiike (1988), and 
Niikura and Nagami (1990). Fast ion distribution functions are used in analytical 
form, which are obtained by solving linearized Fokker-Planck equations. The de- 
pendence of < <ru > on the fast ion energy and ICRF heating is also addressed. 
Bittoni, Cordey and Cox (1980) have shown that the (kinetic) temperature of the 
tail is dependent critically upon the velocity dependence of the particle and energy 
loss mechanisms. Mirin and Tomaschke (1982) have compared the fusion reactivity 
in TFTR, obtained by solving the Fokker-Planck equations for bulk ions with that 
obtained under the assumption of Maxwellian distributions. The effect of the dis- 
tortion of the bulk ion distribution firom the Maxwellian during heating by energetic 
particles was found to be rather insignificant. The heating of plasma by energetic 
particles including loss models have also been studied by Sigmar and Joyce (1971), 
Tsuji, Katsurai, et al. (1976), and Jassby (1977). Steady state calculations were 
made under neutral beam heating scenarios in the presence of velocity dependent 
loss mechanisms (Matsuura, Nakao and Kudo, 1992; Matsuura, Nakao, et at, 1993). 
It is pointed out that some of the calculations are strongly dependent on the plasma 
conditions and the loss mechanisms. 


1.3 Adaptive Node Distribution 


The topic of adaptive grid refinement is wide ranging. The status of the adap- 
tive refinement techniques in finite difference, finite element, and other methods is 
reviewed in, for example, Hawken, Gottlieb and Hansen (1991) and Brebbia and 
Aliabadi (1993). The main steps in the process of grid refinement are: 
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i. the choice of the vaxiable(s) and the error norms based on which the refinement 
is made, 

ii. the computation of error estimates, 

iii. the strategies used to refine the grid. 

These steps are common to both the finite difference and the finite element methods. 
The central idea is to improve the quality of numerical solution by distributing the 
nodes optimally. In this process, the errors are either reduced, or distributed more 
uniformly for the same number of grid points. 

Grid refinement is still new in the area of numerical solution of Fokker-Planck 
equations. The moving finite element method has been applied to solve a single 1-D 
Fokker-Planck equation by Harrison (1988). He has shown that an adaptive grid, 
with less number of nodes, reduces the oscillations in the solution. In the grid re- 
finement process, the error measures are sometimes chosen heuristically, depending 
upon the nature of the problem. For problems in which the variables experience 
steep gradients, the nodes may be placed in proportion to the gradients (for ex- 
ample, Russell and Christiansen, 1978; Dwyer, Kee and Sanders, 1980; Carey and 
Oden, 1984; Amey and Flaherty, 1986). In some other problems, variation in a 
key variable or a parameter relevant to the problem of interest is chosen as the 
error estimate (Yasar and Moses, 1992; Golias and Tsiboukis, 1993; Brebbia and 
Ahabadi, 1993). Altas and Stephenson (1991) have used an approach based on the 
errors in the integral values of the solution. 

Local discretization or truncation errors at different grid points have been used 
for grid refinement in the finite difference method by Crowder and Dalton (1971) 
and de Rivas (1972). Denny and Landis (1972) solved a two point boundary value 
problem adaptively using the leading truncation error term as an error measure in a 
three point finite difference approach. Various strategies for adapting the meshes for 
boundary value problems have been presented in detail by Russell and Christiansen 
(1978). Chong (1978) solved a class of parabolic differential equations using a 
variable mesh finite difference method. The truncation error is distributed more 
or less uniformly. Nonlinear partial differential equations have also been solved by 
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adaptive finite difference method (Schonauer, Raith and Glotz, 1981). IVuncation 
errors have been obtained in a curvilinear coordinate system by Lee and Tsuei 
(1992). Some more details in using the adaptive grid in the finite difference method 
can be found in Berger and Oliger (1984) and Amey and Flaherty (1986). 

The error measures can also be estimated by using a lower and next higher 
order finite difference formulae (Schonauer, Raith and Glotz, 1981). They have 
used Newton’s interpolation polynomial to obtain the coefficients (weights) of the 
difference formulae. Efforts have also been made to generate the finite difference 
weights in a nonuniform grid for derivatives up to any order of accuracy in one 
dimension (Fornberg, 1988, 1990). The error measures can also be obtained by using 
fine and coarse grids. Altas and Stephenson (1991) obtained the error measures by 
dividing the sub cells of the computational domain into smaller cells. In the finite 
element method, the difference between quadratic finite element solution and the 
linear element solution can be evaluated and used as an error measure (Hassan, 
Probert, et ai, 1993). Golias and Tsiboukis (1993) have suggested some more ways 
of estimating the errors. 

Many algorithms have been developed to refine the numerical grids based on the 
error estimates. Chong (1978) presented an algorithm to generate initial nonuniform 
mesh and to refine the grid later by distributing the errors nearly uniformly. Coyle, 
Flaherty and Ludwig (1986) studied the stability of several mesh moving schemes 
that are based on equidistribution. They have studied the schemes that are designed 
to move a mesh so that a particular quantity is uniform over the domain. A simple 
technique has been presented by Sanz-Serna and Christie (1986) to adapt the grid. 
Daripa (1992) developed some grid generators to distribute a fixed number of grid 
points appropriately. Fraga and Morris (1992) also implemented a new method for 
spatial grid refinement. 


1.4 Outline of the Thesis 

The thesis consists of seven chapters. Chapter 2 presents the dimensionless multi- 
velocity scale kinetic equations. The Fokker-Planck collision operators without 
singularities at the boundaries are also given. Chapter 3 deals with the numerical 
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solution of one dimensional Fokker-Planck equations by the Galerkiu finite element 
method and the finite difference method. The computational aspects covered are 
the use of higher order schemes for quadrature and interpolation, conservation of 
density and energy, and the refinement of velocity scales. Redistribution of a fixed 
number of nodes, based on local discretization errors is dealt in Chapter 4. A 
scheme based on the energy error in the computed solutions, to accelerate the time 
marching, is also discussed in this chapter. Solutions of 2-D Fokker-Planck equations 
by the finite element method using five types of quadrilateral elements are given in 
Chapter 5. This chapter also includes the 0-method, and the alternating direction 
implicit method as applied in the finite difference method. Chapter 6 presents some- 
applications involving an improved reaction loss model, and the other sources and 
sinks considered for the simulation of plasmas. Chapter 7 summarizes the results 
obtained in Chapters 2-6 and presents a few suggestions for extending the present 
work further. 



Chapter 2 


Governing Equations and 
Boundary Singularities 


2.1 Introduction 

In a plasma the interactions at the Coulomb collisional time level largely affect 
the energy transfer among the species in order to bring the plasma into kinetic or 
thermodynamic equilibrium. Therefore kinetic models based on the Fokker-Planck 
Coulomb collision operator are essential in representing the plasma behavior. The 
collision operators which model the Coulomb interactions among the charged species 
have interesting properties. They conserve density, total energy and momentum of 
the species. Further the operators vanish for Maxwellian distributions at a com- 
mon temperature in a multi-component plasma. Fokker-Planck equations, when 
represented in phase space, are functions of seven independent variables. In general 
the equations are reduced to 2-D in velocity space under the assumption of spatial 
independence and azimuthal s)Tnmetry in velocity space (uniform magnetic field). 
This is usually done to make the problem tractable. In many applications the 2-D 
models represent the plasma behavior adequately. The equations can further be 
reduced to 1-D in velocity space to deal with isotropic situations. 

When Fokker-Planck equations are applied to model a multi-component plasma, 
the equations are solved for all the general species whose distributions are arbitrary. 
That is a system of coupled nonlinear integro-differential equations is solved. The 


11 




12 


CHAPTER 2. GOVERNING EQUATIONS 


contributions jErom different species are included in the Fokker-Plauc k ccx’ffic ieuts 
through the Rosenbluth potentials and their derivatives. In case of 2-D problems 
the distribution fiinctions, and consequently the Rosenbluth potentials, (,an be ex- 
panded in Legendre polynonaials to obtain the potentials and their derivatives at a 
faster rate. However if some of the distributions are Maxwellian the Fokker-Plauck 
coefficients get simplified. 

In this chapter, the multi-velocity and time scale Fokker-Plauck equations have 
been presented. In this approach dimensionless velocity and time scales rtiost ap- 
propriate for each species are used. The structure of the Fokker-Plauck e(iuations 
in which the coupling between different species is only through the Rosenbluth po- 
tentials makes this approach possible. Analytical limits of the Coulomb collision 
operator at the boundaries are also presented. 


2.2 Basic Equations 


2.2,1 Preliminaries 


In general the distribution functions for charged particles in a plasma are obtained 
by solving Boltzmann transport equation with Fokker-Planck collision teruLS, known 
as Fokker-Planck equations. These equations can be written as 


dt dr nia dv 



(So + Lai 


( 2 . 1 ) 


where /a(r,v,t) is the distribution function for species a at time t and at spatial 
location r with velocity v, F is the force field, (dfjdt), is the Fokker-Planck collision 
term representing the rate of change of due to small angle Coulomb collisions, 
and Sa and La are the source and loss terms respectively. 

Assuming that /„ is independent of spatial position and that force F arises 

only from a magnetic field about which /„ possesses azimuthal symmetry, cq. (2.1) 
reduces to 
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Under these assumptions the Fokker- Planck collision operator derived by Rosen- 
bluth, MacDonald and Judd (1957) is of the form 



1 

2 dvidvj 



(2.3) 


where ga and are called Rosenbluth potentials. The potentials are given by 
(McCoy, Mirin and Killeen, 1981) 

9a = S (^) lnA„(,y/6(v')|v-v'|dv', (2.4) 

(^) fb{y') |v - V\-^ dV, (2.5) 

where the index b is over all the species including a and A^h is the Coulomb logarithm 
which depends on both the interacting species a and b. Coulomb logarithm is taken 
from Futch Jr., Holdren, et al. (1972), which is given by 


In Aab = In 


Triami, \ SxcoacAx) 


,ma + mb 


) 


max 



1 

2 ’ 


( 2 . 6 ) 


where a = 1/137 is the fine structure constant, c is the velocity of light, E is the 
mean energy of the species a or 6 and is the Debye length. Fa and \d are defined 
as 


r„ 


\ 



(2.7) 


where 2’a, rria are the charge number and mass of the species a, rig and Eg are 
average electron density and energy respectively, and eo is the permittivity of free 


space. 
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2.2.2 Basic Equations in 2-D Velocity Space 

Rosenbluth, MacDonald and Judd (1957) further simplified eq. (2.3) by choosing 
spherical poleir coordinates in velocity space (v,d,(p), where v is the speed, 0 the 
pitch angle between velocity vector and magnetic field, and (p is the azimuthal angle. 
With the assumption that the distribution functions / are azimuthally synunetric 
i.e., independent of <p, eq. (2.3) becomes (McCoy, Mirin and Killeen, 1981) 


f^] =r 

[dtj^ dv v^sinO d0 J’ 

where 


if. = + 

The coefficients A^ through Fa are 

2 dv^ ^ dv^ dv ^ dv v dO'^ 

, I ^ 9a _ ^tddQa cots B^Qa 

2 dvdO'^ V ~d9'^ 2 dvd$' 

B^Qa 

2 at;2’ 


Aa = 


Ba = 


C = -1^ + 1-^ 

2vd9'^2dvd9' 

Da = — ^ 4- ^9a sin 9 B^Qa 

2t?2 d9^ 2 dv‘^d9 v 8vd9 

- ^ ^9a , _ . ^dha 

2v^ sin 9 39 2v^ 39^ ^^^^39 

[ 2Dae 2atia9j’ 


( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 

(2.14) 

(2.15) 
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F = sin 6 dga 

The density Ua and energy Ea of species a are defined as 


■a = 27r if fa{v, 9) sin 9 dvd9, 

■aEa = Jo J dvdB. 


(2.16) 


(2.17) 

(2.18) 


In the above equation Ea denotes the average energy per particle. This should not 
be confused with the symbol Ea in eqs. (2.10) and (2.15) which denote one of the 
Fokker-Planck coefficients. The Fokker-Planck equations (2.2) eire now of the form 


^ = r I 1 9Ha \ 

dt “ dv u^sin^ d9 ) 


+ 5a + La, 


(2.19) 


which is a system of equations to be solved for different species in a plasma. 


2.2.3 Rosenbluth Potentials and Legendre Polynomials 


The integrands of the potentials Qa and ha in eqs. (2.4)-(2.5) get simplified under the 
assumption of azimuthal symmetry. The integrands are reduced and expressed in 
terms of complete Elliptic integrals (Rosenbluth, MacDonald and Judd, 1957). The 
potentials are then written in the form of double integrals. As such, in principle, 
these double integrals and their derivatives can be evaluated numerically in order 
to calculate the coefficients. When the potentials are expressed as 




ft. = E 

b 


nia + mb ( 

TUb \Za) 


In Aabhb^ 


( 2 . 20 ) 


( 2 . 21 ) 
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§ 1 , and h(, satisfy the Poisson, equations (Rosenbluth, MacDonald and Judd, 19o7, 
Mirin, 1975), namely 

= - Stt/j,, (2.22) 

= -47r/6. (2-23) 


The Poisson equations may be solved directly to get Qa and ha- The potentials and 
their derivatives can then be obtained from these numerically. 

However faster methods are preferable as these potentials are used repeatedly 
in the calculation of Fokker-Planck coefficients. A method based on which some 
of the Fokker-Planck codes have been developed is to decompose the distribution 
functions and thereby the potentials in Legendre polynomials. This procedure has 
been adopted in the present work also. Writing /j (Mirin, 1975) as 


fb{v, = ^ V^{v, t) Pi(cos 9), (2.24) 

1=0 

where 


Vi^(v, t) = — f 9, t) Pi(cos 9) sin 9 d9. 
l Jo 


(2.25) 


The potentials Qa and ha are written cis 

oo ' 


ga{v,e,t) = lnA„65|’(t;,t)P,(cos0), (2.26) 

ha{v,9,t) = E (^) lii^4(v,t)P/(cos0), (2.27) 


with 


A\(v, t) = 




A-k 

2Z -|- 1 

47r 


r , TOO 

i v^(„\ t) dv' 




, (2.28) 


4/2 
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+ 






I + 1 (t>0 V 


dv' 


(2.29) 


The derivatives of Qa and ha with respect to v can be obtained analytically from 
these coefficients. The 6 derivatives are nothing but the derivatives of Pi{cos6) 
which can be done easily analytically. Legendre polynomials of order up to 8 and 
their derivatives are given in Appendix A. 


2.3 Fokker-Planck Equations in Nondimensional 
Form 

In this section, the normalization constants for velocity, time, and the distribu- 
tion functions are introduced. Making use of these characteristic velocities and 
t im e constants, the Fokker-Planck equations are receist into a nondimensional form. 
The potentials, and density and energy integrals are also written in terms of these 
dimensionless variables. 


Velocity 

In a plasma, whether the species axe Maxwellian or not, ions and electrons have to 
be represented adequately in velocity space for accurate modeling. The cut-off for 
velocity t; is chosen in terms of the thermal speed t;*'* for each species a. Thermal 
speed is defined as 


,.th 



(2.30) 


where Ta is the temperature of the species. In general a minimum of 5 or 6 times the 
thermal speed can provide a satisfactory limit for the representation of species in 
velocity space for many applications. In the present work the velocity normalization 
constant Va is generally chosen to be unless otherwise mentioned. 
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Time 


Here only two time constants, viz. one for all the ion species and the other for 
electrons, are considered based on the following facts. In a typical plasma ions and 
electrons, both fast and slow, have different self relaxation times depending upon 
their energy. The fastest relaxation time is of electrons, tg, compared to ti of ions. 
Among ions, the difference in due to their masses is not significant. Therefore, one 
time constant ii, which is common for all the ions can be effectively used without 
any loss of generality. The self relaxation time U for ions (Roth, 1986) is 


ti 


= 2.256 X 10"‘‘ 


ni(1020) 


s, 


(2.31) 


where Ai is the ion mass in terms of proton masses, the ion density n, is in units of 
10^ m~^, and the kinetic temperature Ti is in units of keV. The time normalization 
constant ii can be chosen as a suitably weighted average of relaxation times U for 
various species. Similarly for electrons the self relaxation time is given by 


te 


= 5.264 X 10~® 


ne(1020) 


s. 


(2.32) 


Here T, is the kinetic temperature of electrons in keV, and n, is the electron density 
in umts of 10^ m“®. In this case tg is taken to be equal to 


Normalizations 

The dimensionless velocity variable and lime v. for each species are defined by 

the following relations: 


V 


VaX 


a ? 


t 


tjTa, 


j = i, e 


(2.33) 

(2.34) 


where and are the normalized velocity and time respectively. The maximmn 
imension ess velocity can be any value, and need not be common for all the 
species in the plasma. It is however convenient to take = l.o for all .species, 



2.3. FOKKER-PLANCK EQUATIONS IN NONDIMENSIONAL FORM 


19 


and choose and change the cut-off velocities during the evolution accordingly. In 
the present work, Xr^ax is taken as 1.0 for all calculations. 

The normalized distribution functions fd^a, 0, Tq) are defined as 




9, T„) = — fdv, 9, t), 

ria 


(2.35) 


where is the density of species a. For simplicity, the same symbol is used for the 
nondimensional distribution function. Now the normaUzed Maxwellian distribution, 
which is independent of 0, is of the form 


faM = 



(2.36) 


where h = mavl/2Ta. The distribution functions for deuterons, tritons and alpha 
particles at a temperature of 20 keV are shown in Fig. 2.1a for a common velocity 
normalization constant u = 6.0 x 10® m s~^. If v is increased further, the spectra of 
all the species become more and more peaked at the origin. This implies that the 
number of particles at the tail (at higher energies) will be negligibly small, which 
can be ignored. The presence of such low values of / at the tail end can become a 
source of oscillation in the numerical solutions. That is, the spectrum may become 
negative at the tail in a numerical solution. 

The undesirable portions in the tails can not be skipped if the same v is used 
for all the species. However, if multi- velocity scales are used, each species including 
electrons can be represented individually, in which case the distribution functions 
of all the species, if Maxwellian, will be as depicted in Fig. 2.1b. It can be seen that 
the number of particles in the high energy region of different species are represented 
adequately and uniformly on their respective velocity scales. 


2.3.1 Fokker-Planck Equations 

The new variables incorporated into eqs. (2.8)-(2.16) lead to the following multi- 
velocity and time scale kinetic equations in dimensionless form. Same S 3 Tnbols have 
been used for the normalized quantities. This has been done for simplicity and to 
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a) Common velocity scale 



Figure 2.1; Maxwellian distributions for different species 
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maintain continuity. Henceforth, unless otherwise specified, all the references will 
be made to these dimensionless variables only. The Fokker-Planck equations in 
nondimensional form are 


dTj “ x^sin^ 86 J 


■j- Sa -h La, 


where 


Ka = 




j = i,e 


and 


Ga 

Ha 


dxa 86 


Dafa + Ea 


dXa 


+ Fa 


86' 


(2.37) 


(2.38) 


(2.39) 

(2.40) 


Sa and La are the source and loss terms, respectively. The Fokker-Planck coefficients 
are 


^ 1 ^9a 

^ 2 8x\ 5x2 ^°8Xa Xa 86^ 

1 d^Qa _ cot 6 dga COt^ 8'^ga 

2 8x^86^ Xa 86 2 8xjd6' 

\ - 

“ 2 8x1' 

' == 1 1 ^9a 

“ 2xo 86 2 8xJ86' 

^ _ sivi68^ga sin 6 d^ga sin0 8^ga 

“ “ 2 5x25^"^ Xa 8x^86 

1 dga ^ cos 6 8^ ga . ^^ha 

2x2 gJjj 0 Q0 2x2 Q02 Q0 ’ 


(2.41) 

(2.42) 

(2.43) 


(2.44) 
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1 dga 1 

2x^ de 2dx^d9 


sin 9 d^Qa sin 9 dgq 
2x1 89'^ 2Xa dxa 


(2.45) 

(2.46) 


The Rosenbluth potentials g^ and ha, S'nd their derivatives are now of the form 


ga{Xa, d,T) 


(2.47) 

dga 

dx^ 

= §?(!)“ 

(2.48) 

d^ga 

dxl 


(2.49) 

d^ga 

= EE(|)^nA.(^)(|)^f .„co.., 

(2.50) 

ha{Xa, 9, r) 

;^Y mj \Za) ‘^\naJ\vJ 



A'l{Xa,T) Pi{cOs9), 

(2.51) 

dha 

dx^ 

■y-' y-v "b '^b f f f 




(2.52) 


The expressions of A^(xa, r) and Bi{xa, r), can be written in terms of four function- 
als Ml, Nj, Ri, and Ei, as follows: 

k 4'7r 

(2.53) 

47r f 1 


Bl(Xa,T) = 


2Z -h 1 ( (2/ -h 3) 
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where 


(2.54) 


Va 

U = —Xa. 
Vb 


(2.55) 


VJ*”s are as given in eq. (2.25), and the functionals are defined as 


Ml (W (Xa)) = 

f'^w(y)y^^-Uy, 

Ju 

(2.56) 

Nl (w (Xa)) = 

f w{y) dy, 

Jo 

(2.57) 

Rl (W (Xa)) = 

nw{y)y^^-Uy, 

Ju 

(2.58) 

El (W (Xa)) = 

f w(y)y^*'^^^ dy, 

Jo 

(2.59) 

with 



_ ^ 

^max — ^max- 

Vb 


(2.60) 

It may be noted that maximum limit of u (tWc) 

can be taken as Xmax itself. Since 

Xmax has been chosen as 1.0 in the present work, 

Umax also Can be taken 2is 1.0. 


As mentioned earlier, the derivatives of Aj and Bf in terms of the new variables 
can be obtained analytically. They axe 


dA\ 

dx^ 

dxl 


= I u'-'M. {Vh - (! + 1) u-’-^N, (V;‘) ], 

= (I + m + 2) (V;"') + i (; - 1) J-=M, (v;‘) ) 


( 2 . 61 ) 



- Aw V*{u), 

= 2i^ { (21^1 > 


(2.62) 
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d^Bf 


dxl 



dxi 


(V1‘) - (i - 1) u-‘N, {V,‘) ] 

[21 Ij 


Atv 

21 + 1 


t (2Z + 3) 


lu-‘'^Ei(V,'’)+v‘Mi(V,‘’)] 


(V,") + (v;‘) i| , 


(2.63) 


(2.64) 


- 2) u‘->je, (V^‘) -(l + l) u-‘-^N, (v;‘) ]| , (2.65) 


+(/ + 1)(Z + 2)(Z + 3)(Z + 4) u-^-^Ei (VJ') ] 

" (2r^f ^ ~ 

+(Z - 1) Z (Z + 1)(Z + 2) u-‘-^Ni (Vi’’) ]} 

-Stt VJ^(tt). (2.66) 


With the normalized distribution functions, the density and energy integrals in 
eqs. (2.17) and (2.18) become 

2. rr /,(...«) aJ^sin^tZxad^ = 1, (2.67) 

Jo Jo ®aSin6' cZaj^cZ^ = (2.68) 


where Ea is the energy per particle. 
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Boundary Conditions 

The Fokker-Planck equations (2.37)-(2.46) are solved with the following velocity 
space boundary conditions (Killeen, Kerbel, et al., 1986). 


faixa = 0, 9) is independent of 6, 

(2.69) 

^(a;, = O,0 = 7r/2) = O, 

(2.70) 

QQ {Xa, 9 = 0)= -^{Xa, 9 = T^)=0, 

(2.71) 

fai^a ” ^max') ~ 0. 

(2.72) 


Initial Conditions 

Initial distribution of each species fa(Xa, 0) at time Ta = 0 is specified. can be any 
distribution including non-Maxwellian. This would facilitate a study of relaxation 
of arbitrary distributions to approximate steady state Maxwellian in the absence 
of source and loss terms. In fact this can serve as a test case for the code(s). In 
general, the initial distributions are used in a form given by (Futch Jr., Holdren, 
et al., 1972; McCoy, Mirin and Killeen, 1981) 

fa(Xa, 0) = eta exp {-ba{Xa - f3af “ Ca(0 - 4)^), (2-73) 

where 4 is defined as = IrUaPlIvl and aa,ba,Ca, and 4 are constants which 
are given as input. is obtained by satisfying eq. (2.67). The initial distribution 
functions of the form given here are referred to as Gaussian distribution functions. 
When Ca = 0, i.e. /„ is isotropic, the above distributions are sometimes called as 
shifted or displaced Maxwellians. Further when the distributions are Maxwellian, 
some of the Fokker-Planck coefficients also get simplified. 
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2.3.2 Background Species 

Siinulatioiis of plasmas include both energetic and warm/cold ba( kgroiind popula- 
tions. Energy from fast ions is transferred to heat the background sp(K:ies. The 
background species are in general assumed to be isotropic, that is there is no 6 
dependence, and also Maxwellian. The distributions are of the form 


fix) = 



(2.74) 


where b = mv^/2T and v is the velocity normalization constant. In this case, the 
potentials Qa and ha get reduced to simpler forms and the contribution from such 
distributions to the respective coefficients are 


Bt 


mi, 


(^y) InAofc — 


\ZaJ ria 

V TT 2 


mb Vb 1 ,b 
ma Va 2lm 


ZbY, . rib sinO 
InAoi — - — 
ria 2%a 


(I) 


1 - 




<i){\/hu) -t- 


sfr^Xa 




(2.75) 

(2.76) 

(2.77) 


where u is as defined in eq. (2.55), with Vb = v. In eqs. (2.75)-(2.77), the subscript 
or superscript b refers to the background species under consideration. In the above 
expressions terms containing e have been ignored because they are negligible if 

suitable cut-oflF velocities are chosen. The function <j>{z) is the error function defined 
as 



(2.78) 


which can be evaluated easily by power series approximation. For z < 0.6, 


k=0 


^2k+l 

k\{2k + l)' 


(2.79) 
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and for z > 0.6, 
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<l){z) = 1 - 




(2.80) 


with 


77 = 


1 

1 +pz' 


(2.81) 


and p = 0.3275911, Ai = 0.225836846, A^ = -0.252128668, Ag = 1.25969513, 
A 4 = -1.287822453, and Ag = 0.94664607. 


2.4 Fokker-Planck Equations in One Dimension 


One dimensional Fokker-Plandc equations are easily obtained by setting Z = 0 in 
eqs. (2.37)-(2.46) and neglecting terms containing 9 derivatives. This results in the 
following equations: 


dfa ^ 1 dGa 

dTa "Xldx^' 

^ = 7^1 A 

dTa "xldXa 

The boundary conditions are 

= 0 at Xa = 0, 
ffi 0 at Xfl X m/iT — 1.0. 


Aq + Bf 


dh 

dXa 


(2.82) 

(2.83) 


(2.84) 

(2.85) 


The coefficients and the potentials are 


A. 


2 dxl dx^ dXa °'dXa 


( 2 . 86 ) 
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(2.87) 

(2.88) 

(2.89) 


The expressions for and Bq are given by 

A>(x^,t) = 47r[i«-o(/) + *fo(/‘)l, (2-90) 

L U 

SJ(x„t) = 47r[i(ia(/‘)+“'M.(/‘))+(«.iVo(/‘)+ft(/‘))]- (2.91) 


The derivatives of Ag and Bq are 


dAi 

dx^ 

= [-!«,(/*)], 

(2.92) 

8x1 

= 4x f (/‘)1 - 4,rA(u), 

(2.93) 

a-Bj 

dx^ 

= 4x [ i (2uAf„ (/‘) - ~Eo (/*)) -t- No (/‘)] , 

(2,94) 

d^Bl 

8x1 

= ^’^[|(;^^»(rt + "»(/*))], 

(2,95) 

8=BS 

dxl 

= 4x[-1eo(/‘)], 

(2.96) 

d^Bl 

dxl 

= 4’^[fS)(/‘)]-8xA(u). 

(2.97) 
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The functionals Mo,No, and Eq are defined in eqs. (2.66)-(2.59). The coefficients 
can further be reduced to 


“ \ZaJ KriaJ rrib 

= 4' E (#)' l-A.* (^) (?=)' 

\naJ \VbJ 




The density and energy integrals in 1-D are given as 

rl 

AtT I fa{^a) — 1 , 

Jo 

'^aVa f ^ 0 (^ 0 ) dx^ — E^, 

Jo 


(2.98) 

(2.99) 


( 2 . 100 ) 

( 2 . 101 ) 


where Ea is the energy per particle, as earlier. As in 2-D, if the distribution is 
assumed to be Maxwellian for a background species, the corresponding contributions 
to the coefficients can be simplified. The final form of the coefficients are the same 
as given in eqs. (2.75)-(2.76). 


2.5 Collision Operators without Boundary Sin- 
gularities 

The Fokker-Planck equations both in 1-D and 2-D can be solved numerically by the 
finite difference and the finite element methods. In the ceise of the finite difference 
method, the singularities in the collision operators at the boundaries have to be 
dealt with properly. The usual method employed is to approximate the operator 
at a point close to the boundary. In some ceises, especially when the ntimber of 
grid points in the computational domain is large, the boundary point can as well 
be ignored and the solution at the boundary is set to the value at the nearest node. 
In this section the analytical limits of the collision operator at the boundaries are 
presented, which are used in later chapters for discretization at the boimdaries. 
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2.5.1 1-D Operator 

The 1-D Fokker-Planck equations (2.82)-(2.83), dropping th('. subscript a and omit- 
ting the constant Ka for simplicity, can be written as 


df _ 

dr dx ’ 

Af + B^ 


dr x^ dx 


( 2 . 102 ) 

(2.103) 


Here G is zero at a; = 0 because A = B = 0sLtx = 0. It can be seen oiisily that 
dG/dx and d^G/dx'^ are also zero at x = 0, as must be the case if the limit of the 
right hand side is to exist. Therefore the limiting value can be found by applying 
L’Hopital’s rule (Apostol, 1962). Before doing this it is convenient to rewrite': the 
right hand side of the above equation in expanded form as 




(2.104) 


where the coefficients A\B*, and C*, obtained from eqs. (2.86)-(2.87) are given 
below: 


2 dx^ X dx^ X dx dx^ ’ 
dx^ x dx^ x^ dx dx ’ 
2dx^' 


It should be noted that at x = 

a = ^ = ^ = h-^ 
dx dx^ dx 


0 , 

= 0 . 


(2.105) 

(2.106) 
(2.107) 


(2.108) 


Therefore, second and third terms in A* are of the form 0/0. Their values arc 
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lim 

x— ♦O 


X dx 


dx^ 


Thus, the coefficient A* at x = 0 is given by 


( 2 . 110 ) 


2 dx^ dx'^ ' 


( 2 . 111 ) 


It is easily seen that S* — > oo, as x — »• 0. This is due to the third term in eq. (2.106). 
Thus the limit of the second term in eq. (2.104) must be obtained together with the 
boundary conditions (2.84). Its limit is given by 


Um B 

x — »0 


dx 


dx^ dx^ 


( 2 . 112 ) 


This term is then combined with the third term in eq. (2.104) which is finite and 
nonzero at x = 0. The resulting operator at x = 0 can be written as 


dr 


= A*f + C 


dx^ 


(2.113) 


where 

2dx^' 

and A* is as obtained in eq. (2.111). 


(2.114) 


2.5.2 2-D Operator 

In this case the operators have singularities at three boundaries, x = 0, 0 = 0, and 
9 = TT. The singularities at ^ = 0 and tt are of first order, due to the division by 
sin0 term, while the singularity at x = 0 is of second order (see eqs. (2.37)-(2.46)). 
As the coefficients Aa through Fa are evaluated using potentials Qa and ha which 
involve Legendre polynoimals Jf^j(cos0), the analytical limits must be valid for all 
the orders, Z, of Legendre polynomials. 

The general steps are similar to the case for 1-D Fokker-Planck operator, but 
considerably more involved due to large number of terms and the need to evaluate 
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the limits for all values of I. Analysis of singularities is somewhat fadlitato<i hy 
rewriting eq. (2.37) in the following manner; 

df _ 1 ^ 1 djsmeH*) (2.115) 

dx^ x^sinO 66 


where 


G 



Af + B%^C 


69' 


D*f + E*^ + F 

OX 


66 ’ 


The coefficients D\E% and F* are given by 


(2.116) 

(2.117) 


I 6^g I 1 6i^9 ^ ^9 

^ ~ ^66^'^ 2 6x^66^ X 6x66 2x'^sm^9d6 

cot 6 6^g 6h 

'^'2^W~d6' 

2x69 2 6x66' 

F* = + 

2x^ 69'^ 2x 6x 


(2.118) 

(2.119) 

( 2 . 120 ) 


In the above equations, the constant Ka and the subscript a have been omitted 
for simplicity as earlier in this section. The following points have been mad('. ilsc 
of in the process of obtaining analytical expressions for the 2-D operator at the 
boundaries; 


a. At X = 0, coefficients A-C, D* — F* are zero. 


b. At X = 0, 5 ’ and h and their derivatives with respect to x are zero for I > 4. 

c. At X = 0, the terms containing 6 derivatives are zero. 


d. At 0 = 0 and tt, A ^ 0, B ^0, and C = D* ~ E* = 0. 
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e. All the odd derivatives of Pi{cos9) are zero at 0 = 0 and tt. 

f. The boundary conditions df /dO = 0 at 9 = 0 and tt axe used where necessary. 


At 0 = 0 and tt 


Considering the above points, the term containing df/d9 in G drops out at 0 = 0 
and TT, while the other two terms remain, with the l imi ting expressions for A and 
B as follows: 


A = d^g 

2 dx^ dx^ dx ^ dx^ xde^'^dxdB^' 

B = — ^ 

2 5x2- 


( 2 . 121 ) 

( 2 . 122 ) 


The limit of the second term in eq. (2.115) is same at 0 = 0 and tt, and is given by 


2dH* 2 fdD* dE*df 

x2 89 “ x2 QQ •/ + Qg Qg2 J ’ 


(2.123) 


where the terms which axe zero have been dropped. The expression F* remains the 
same as in eq. (2.120). The limiting expressions for the remaining two coefl&cients 
are as given below: 


dD* 

2 d^g 

1 d‘^g 1 ^g 1 d^g 

02/l 

89 


3x2 g. Qxd9'^ 2 0x2002 

002’ 

dE* 

1 

, 1 ^0 


89 

2x002 

"^2 0X002' 



(2.124) 

(2.125) 


At X = 0 


Equation (2.37) is expanded and written in the following form to facilitate the 
analysis of the singularity at x = 0: 


dr 


.df -d^f - df 


+ 


(2.126) 
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where 

A 

B 

C 

D 

E 

F 




dA 1 dD 
H — r 


dx sind dO 


x^ 

1 

x^ 


\ dB 1 dE 
A+ — + - 


dx sia0 dO 


B, 


x^ L 


C + 


sin0 


E 


x‘- 


ax sin0 sin 0 50 


x^ Lsin0 J 


( 2 . 127 ) 

( 2 . 128 ) 

( 2 . 129 ) 

( 2 . 130 ) 

( 2 . 131 ) 

( 2 . 132 ) 


Making use of the observations made earlier, it can be easily seen that the tc'rnis 
contaimng B, E, and F drop out at x = 0, because their limits are z<'r<). As in 
1-D, B -> cx), as X 0. However, if df/dx is taken as zero at x = 0 for all values 
of 0, the limit of the term containing B exists, and can be merged with th(‘ term 

containing C, which is finite at x = 0. The resulting Fokker-Planck op<^rator at 
a; = 0 is: 


dr + 


( 2 . 133 ) 


where .4- and C' are the same as in 1-D, eqs. (2.111) and (2.114), It i., .seen 
that the Wing expressions of the operator in 2-D is same as in 1-D. It nnrst be 

mentioned here that for the limit to exist at x - 0 in 9 n fK k i 
(a n L at X - u in 2-D, the boundary condition 

I ';” : " of at 9 = ,r/2 only 

norasf ,r: ~ - --‘■y 

o^rTiitr 

Ot tne collision operator at the boundaries. 
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2.6 Summary 

The work that has been carried out in this chapter can be summarized as under: 

i. 2-D velocity space kinetic equations applicable to plasmas in their general 
form have been presented. 

ii. The concept of using multi-velocity and time scales in the kinetic simulation 
of plasmas has been presented. The dimensionless Fokker-Planck equations 
and Rosenbluth potentials in 2-D velocity space are given. These equations 
can be used for an arbitrary number of general and background species. In 
the Ccisq of background species, separate expressions for the coefficients have 
been arrived at using the characteristic velocity scales. 

iii. One dimensional multi-velocity scale Fokker-Planck equations, as derived from 
2-D equations, are also presented. 

iv. Both 1-D and 2-D collision operators have singuleirities at the boundaries. 
The singularities are removed analytically and the non-singular operators to 
be applied at the boundaries are given. 



Chapter 3 


One Dimensional Fokker- Planck 

Equations 


3.1 Introduction 

In a numerical solution of the Fokker-Planck equations, accuracy of the computed 
solution with regard to the distribution functions as well as satisfying the con- 
servation requirements are equally important. Difficulties arise as the numerical 
solutions are expected to conserve both the density and energy. Suitable choice of 
boundary conditions can sometimes be applied to satisfy either the density or the 
energy condition accurately. In this context the finite element method can possibly 
play a better role because it is known, despite the fact that more computational 
effort is involved, that the spatial resolution in the finite element solutions is better 
than the finite difference solutions. In addition, the finite element method offers 
certain advantages for problems with ciuved boundaries. The operator singularities 
encountered in curvilinezir coordinate systems at the boundaries aire also easier to 
deal with in the finite element method. 

In this chapter one dimensional Fokker-Planck equations have been solved by 
both the Galerkin finite element method and the finite difference method. Vari- 
ous aspects of the computation process are illustrated for a multispecies plasma 
including electrons as a general species. The Fokker-Planck solutions have also 
been compared with the Spitzer energy relaxation model which is valid under the 
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Maxwellian conditions. 


3.2 Finite Element Method 

In this method, the computational domain of interest is dividcnl into a uiuuher 
of non-overlapping subdomains or finite elements. The solution to th<‘ governing 
equation is then constructed by using piece wise polynomials over each eUmumt, 
called the shape functions. This approximate solution in satisfying th<' governing 
equations results in a residue R which can be minimized with n'spoct to weight ing 
functions Wj as 



where N is the number of nodes in the domain. In the finite clement (iis<:r<'tizatiou 
procedure, the elemental contributions are assembled to obtain the global system of 
equations at all the nodes. The global equations are solved to get the approximate 
solution to the governing equation. Nevertheless there are many ways to choose the 
weighting functions Wj. In Galerkin finite element approach, the Wj are choscm to 
be identical to the shape functions themselves. 

The various steps associated with Galerkin finite element method of solution are 
given in detail with the help of 1-D Fokker-Planck finite element formulation. 


3.2.1 Galerkin Formulation for 1-D Equations 

The 1-D Fokker-Planck equation to be solved is chosen in its asual form as in 
eq.(2.83) with some additional terms. The additional terms include a term of the 
form Qf and a constant term S. Omitting the subscript a and the constant K, the 
equation is written as 
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• • 

Linear 


• • • 

Quadratic 


• • • • 

Cubic 


Figure 3.1: One dimensional finite elements 


with the boundary conditions 

(3.3) 

(3.4) 


= 0 at a: = 0, 
ax 

f = 0 at X = 1.0. 


Here the computational domain is restricted between 0 and 1.0. Let / be the 
approximate solution to this equation. The residue R can then be written as 


R 


J_A 

dr x^ dx 


Af + B 


dx 


-Qf-S, 


(3.5) 


and its integral with weighting functions Wi is minimized over the entire domain as 


[^'^WRA7rx'^dx = 0, i = l,...,N, (3.6) 

Jo 

The factor 47 rx^ arises out of the governing equation being written in spherical 
geometry. Wi are chosen as follows: The domain is divided into Nx elements, not 
necessarily having equal lengths. The number of unknowns at each node is only one. 
Therefore, linear elements having two nodes are represented by two linear shape 
functions. Similarly higher order elements namely quadratic, cubic etc. have 3, 4 
and more nodes respectively. These elements are represented by the corresponding 
number of shape functions. In the present work, the nodes are placed in each 
element at equal intervals for quadratic and cubic elements as shown in Fig. 3.1. 
The approximate solution f at time r can be written over each element as 

m 

;=i 


(3.7) 
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where m is, the number of shape functions or, the number of nodes in eaeii eleiiKnit. 
Here fj are functions of time and the shape functions dcpemd only on x. Ciioos- 
ing <l>, Is weighting functions, the residual equations for an ek'iruuit o with global 
coordinates and a?e-i-i can be written as 




dr dx 


Af + B 


dp 

dx , 


Qfx^ — Sx 


(l>i{x) dx — 0, 
i = 1, . . . , m , 


(3.8) 


for all the shape timctions m in the element. Integrating the second term by parts 



The contribution of the last term in the left hand side of the above equation vanishes 
at the boundary nodes. That is, it can be seen that 




at a: = 0 and 1.0 , 


(3.10) 


because A and B are zero at x — 0, and / = 0 and H ~ 0 at x = 1. Substitution of 
eq.(3.7) for / in eq.(3.9) leads to 


m 


E 


<f>i (j>j x^ dx 



A<j)j B 


d4>j 

dx 



fAr) 


- (/ M 4>j dx'j fj{T) = x^ dx 


+ ye+l - 2/e, 


i = l,...,m , 


(3.11) 
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with 

(3.12) 

(3.13) 

t^ + uP(T) = v(T), (3,14) 

where 


iv = 

j (pi ipj dx. 


(3.15) 

tiij = 

r^B^dx+, 

J dx dx J 

r r ■) 

1 A <j>j ~ J 4>iQ 

(3.16) 

Vi = 

f <f>i S d/x + — 

Ve+l- 

(3.17) 


The choices of shape functions and elements, and thereafter the computations of 
the element matrices are discussed in Section 3.2.2. The element equations are 
generated for all the N^ elements. Finally these equations are assembled to arrive 
at the global equations. The global equations in matrix form are 

T-^ + Uf=V, (3.18) 

dr 

representing a system of coupled ordinary differential equations. The discretization 
of spatial part of eq. (3.2) by finite element method has resulted in a system of 
ordinary differential equations in time. These equations are integrated in time by 
the one step 0-method. The boundary conditions are implemented as follows: 

1. The derivative condition at a; = 0 is homogeneous and as such does not 
contribute to eq. (3.11). 


Ve 


y-a*. 


x—x^ 


y'» = 


X—Xc+l 


These element equations in matrix form are 
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2. The other condition that / = 0 at a: = 1.0, is imposed before solving the 
global equations. 

3.2.2 Computation of Element Matrices 

Each entry of the matrices t,u and v in eq. (3.14) requires nunu'rical evaluation 
of 1-D integrals. A local coordinate ^ between -1 and +1 is used to transform the 
integrals in order to facilitate the use of Gauss quadrature. Introducing thc^ local 
variable ^ with the transformation 

+ -1<C<+1, (3.1!)) 


with 


Xh •C'e+1 ®ej 


(3.20) 


the equations (3.15)-(3.17) become 

tii = Y (3.21) 

Xh 

“ y /-I (3.22) 

"" y Li + ye - yef 1 . (3.23) 

The 1-D shape functions in terms of the local variable ^ for different elemcmts are 
given m Appendix B. The integrals are evaluated numerically using Gauss cpiadra- 
ture, which, for a general integral is given by 

J ^iO(^ = Y>Wkz{^k), 

k=l 


(3.24) 



3.2. FINITE ELEMENT METHOD 


43 


where and w/t are the Gauss points and weights respectively. The entries in the 
element matrices given in eqs. (3.21)-(3.23) are calculated by the above procedure. 

To evaluate the element matrices, the Fokker-Planck coefficients A and B are 
needed at each Gauss point These coefficients are not in explicit form for direct 
computations. Calculating the coefficients eiccurately at each Gauss point for indi- 
vidual elements at every time step is quite costly and it should possibly be avoided. 
However, the coefficients can be calculated after every time step at all the nodes in 
the domain. Then these nodal values can be used to interpolate for any value of 
within each element. In other words, the coefficients A and B are written simply as 


m 



= 

j-i 

(3.25) 

B’iO 

m 

(3.26) 


in ternas of the nodal values Aj and Bj and the shape functions <f>j. The interpolated 
values axe used in the integrals of element matrices. 

Sinailarly, in the case of source and loss terms, where expUcit expressions are not 
available, e.g. reaction loss, the same procedure has been employed. For example, 
Q can be written as 

m 

= P-27) 

j=i 

where Qj are the nodal values. In the present case, the above procedure is supported 
by the fact that the Fokker-Planck coefficients behave smoothly. 

With the help of these representations, it is now possible to calculate the ele- 
ment matrices and finally the global matrices T, U,V. This completes the spatial 
discretization procedure. Next step is to integrate the equations to advance in time. 

3.2.3 Time Integration 

In the Fokker-Planck equations nonlinearity comes through the potentials which 
are dependent on the distribution functions /. There are different ways to solve a 
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system of nonlinear equations. For example 

i. The coefficients at the n+lth step can be obtained by extrapolating th(' values 
at steps n and n - 1. The resulting linear system is solved. 

ii. The system is treated as quasi-linear. The coefficients are calculated using 
the distribution functions at the previous time step. 


In both the cases, the coefficients can also be iterated over the solutions. In the 
present work, the method (ii) has been adopted because past research(‘rs have found 
that the coefficients of Fokker-Planck equations are slowly varying (Kill<'<ni and 
Futch, 1968; McCoy, Mirin and Killeen, 1981) and in general it is not u<'c<'ssary 
to iterate the coefficients over the solutions. Further, the coefficients lu'e relativtdy 
insensitive to the distribution functions because they involve only their intf'grals. 
Provisions have, however, been made in the computer codes to iterate the co<’fficients 
also. In general, iterations have not shown any appreciable difference in th<‘ results 
obtained. 

The temporal part is discretized by standard finite difference te(Iuu<iiu\s. Th(^ 
scheme that has been employed here is generally known as the one step 0 nw'thod. 
This method is chosen because it gives a chance to conduct many nuim'rical (‘xper- 
iments by varying 0 values. This can also help to study nmnerically th<! stability 
and convergence aspects of the computed solutions. 

The 0 method approximates the time derivative as follows. A w(dghted av<'rag (5 
of the time derivatives of the dependent variable f at two consecutive time st<'ps, 

say n and n-j-l is used to interpolate the values of jf linearly between the two steps. 
That is 


Using this approximation, eq. (3.18) becomes 


jpl4-l 

= e(V“« - ur«) + (1 _ e) (V” - u") , 


(3.28) 


At 


(3.29) 
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or, 

+ ©uj - (1 - ©)u) r + (©v"+^ + (1 - e)V") .(3.30) 

The right hand side can be calculated using / at time step n. Now the resulting 
algebraic system of equations can be written as 

= B, (3.31) 


where 

A = (^T+eu), 

B = - (1 - 0)uj r + (©v"+' + (1 - e)v") . 

Different values of the parameter © lead to the following schemes: 

© = 0.0 forward or explicit method (0 (At)) , 

= 1/2 Cranck-Nicolson method (C>(At^)) , 

= 2/3 Galerkin method (0( At)) , 

= 1.0 backward or fully implicit method (0 (At)) , 
where 0 (At) is the order of convergence of the respective methods. 


(3.32) 

(3.33) 


3.3 Computation of Fokker-Planck Coefficients 

The coefficients which appear in the Fokker-Planck collision operators are calculated 
numerically at the grid points. The collision operators can be written as a sum 



(3.34) 
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where {dfa/dtf, represents collision of species a with species b. Thv. contributions 
from all the general species b including a are accounted fis outlined in ('(js. (2.88) 
(2.97) and that of background species are calculated as given in S<'ction 2.3.2. In 
multi-velocity scale approach numerical approximations that aris(' primarily are 
through the Rosenbluth potentials and their derivatives. The computation.s involve: 

i. numerical evaluation of the integrals M,N,R, and E in oxis. (2. .>6) (2.a9), 
and 

ii. numerical interpolations which arise out of repr<'s<'ntiug sped('s on different 
dimensionless velocity scales and grids. 

The values of the functionals obtained at grid points are. used to calculate th('. 
Rosenbluth potentials and their derivatives. 


3.3.1 Numerical Quadrature and Interpolation 

Various schemes can be used for numerical integration and interpolation. Higher 
order quadrature and interpolation will naturally lead to improv«'nient in the fxccu- 
racy of the coefficients. The limits of the integrals in eqs. (2.56) ■ (2.59) can vixry by 
a factor of 40-60 due to difference in the velocity cut-offs. Since (maximum 
cut-off in the normalized velocity) is always 1.0, all those integrals iu-e (evaluated 
first for each species b using its own uj,. Then the integral values at the grid points 
are interpolated for intermediate x values as required by the spcci(^s a. In this work, 
two schemes have been implemented for quadrature and interpolation. First one 
is the trapezoidal rule which is commonly used and easy to impleimuit. In order 
to improve the accuracy in the calculation of coefficients, a scheme based on cubic 

polynonaials has also been used. In the trapezoidal rule the integral is aj)proximated 
as 


y(®) 


(3.35) 
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Table 3.1: Fokker-Planck coefficient Aa obtained by using different integration 
schemes (S)- E-Error function, T-Trapezoidal and H-Hermite 


Ax 

S 

X2 




xe 

xr 

Xs 

Xg 

1/30 

E 

0.0446 

0.3443 

1.095 

2.393 

4.218 

6.449 

8.901 

11.370 


T 

0.0663 

0.3836 

1.142 

2.437 

4.251 

6.466 

8.900 

11.350 


H 

0.0445 

0.3444 

1.095 

2.392 

4.216 

6.447 

8.898 

11.370 

1/60 

E 

0.0056 

0.0446 

0.148 

0.344 

0.654 

1.095 

1.674 

2.393 


T 

0.0084 

0.0501 

0.156 

0.354 

0.665 

1.107 

1.686 

2.404 


H 

0.0056 

0.0446 

0.148 

0.344 

0.654 

1.095 

1.674 

2.392 

1/90 

E 

0.0017 

0.0133 

0.045 

0.105 

0.202 

0.344 

0.537 

0.786 


T 

0.0025 

0.0149 

0.047 

0.108 

0.206 

0.348 

0.542 

0.791 


H 

0.0017 

0.0133 

0.045 

0.105 

0.202 

0.344 

0.537 

0.786 


In the case of cubic approximation both the Hermite and Lagrangian polynomials 
have been used, with almost identical results. For the interpolation also, both lineeir 
and cubic polynomials have been implemented. It should be mentioned here that 
1-D Fokker-Planck coefficients axe more easily obtained through eqs. (2.98)— (2.99) 
than the procedure outlined in eqs. (2.86)-(2.97). In the case of 2-D coefficients the 
steps discussed above are in principle the same and these steps are repeated for all 
the orders I of the Legendre polynomials. 

3.3.2 Comparison of Accuracy 

In order to evaluate the merits of the trapezoidal and Hermite quawlratures, Fokker- 
Planck coefficients are calculated for a single species with Maixwellian distribution. 
The values of 1-D coefficients A^ and Ba are also obtained accurately by using error 
function given in Section 2.3.2. The coefficients Aa computed by employing different 
schemes are shown in Table 3.1. The values are given for the first 9 grid points 
(excluding a; = 0, where A^ = Ba = 0) corresponding to different mesh sizes Ax. 
It can be seen that Hermite scheme is more accurate even with less number of grid 
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points than the trapezoidal rule. Such high accuracy is maintained in the case of 
multispecies also. Coefficient reaches asymptotic value beyond certain x. So 
both the schemes produce the same results. 

Coefficients Ba which involve the functionals -M and E have not shown any 
significant difference in the low velocity region (near x 0) for both tlu^ s(.h('incs. 
This is because the functional M dominates near x = 0. Con.sequently the int('gral 
varies slowly in the low velocity region. 


3.4 Convergence 

Among the various time stepping schemes presented earlier (for diffen'nt 0 values), 
Crank-Nicolson and backward (fully implicit) schemes are unconditionally stabl(> 
for any value of the time step. However the former, although s(’cond order accurat<^, 
is prone to oscillation in the solutions. Backward and Galerkin sclu'im's are first 
order accurate and the later is only conditionally stable. Therefore the t<'.mporal 
behavior of these time stepping schemes has to be assessed. Similarly tlwi spatial 
convergence behavior of the various finite elements must be studied in order to have 
stable and accurate numericeil solutions. 


3.4.1 Error Norms and Test Input 

An essential part of any computational method is the assessment of the accuracy 
of the results obtained. For problems related to transport phenonu'na, in which 
the conservation of certain quantities is important, density and energy norms of the 
error in the computed solutions are quite often used. The relative error norms can 
be defined as follows; 


Ln norm of A/ 
Ln norm of / 


7nlA/|”dn 


L 


l/n 


with 


(3.36) 


|A/I = l/.„ - 


(3.37) 



3.4. CONVERGENCE 


49 


where /„ and fapp are the exact and approximate solutions respectively. Here four 
types of relative errors namely, 600,61,62, and Crms corresponding to the Loci'!, L2, 
and Lrms norms are defined. For 1-D problems, the errors in percentage are mea- 
sured as 


600 


ei 


62 


6rms 


mi 


fn 


X 100, 


/I A/I xUx 

J fx'^dx 


X 100, 


/ |A/|^ x^dx 
J p x^ dx 


1/2 


xlOO, 


N. 


EiA/.r 


1/2 


X 100, 


(3.38) 

(3.39) 

(3.40) 

(3.41) 


where the subscript i is over all the nodes in the velocity grid. The Ln norms can 
also be used to measure the degree of relaxation of of an arbitrary distribution at 
different times. In this case the represents the Maxwellian distribution to which 
the initial distribution function is expected to evolve. When used for this purpose 
the norms are denoted by L„, and not by e„. 


Test Input 

Numerical experiments have been conducted first for a single species in the absence 
of source and loss terms with Maxwellian and Gaussian distributions. An arbi- 
trary initial distribution (at r = 0) can be expected to relax to an approximate 
Maxwellian in about 6t (t -relaxation time). The following aspects of evolution are 
covered in the numerical tests conducted: 

i. If a Maxwellian itself is given as input, the distribution function should not 
evolve (change) under the action of Fokker-Planck collision operator. How- 
ever the munerical approximation of the collision operator induces numerical 
evolution (Sued, Appert, et al., 1986), and the input steady state distribution 
begins to deviate fi:om the Maxwellian. This occurs due to the contribu- 
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tions from both spatial and temporal discretization proceduros. Keeping the 
time step constant, by varying the number of grid points, the effect of .spatial 
discretization errors on the artiHcial (spuriou.s) numerical evolution can l>e 

studied. 

u. When arbitrary initial distributions (used mostly in the form ofCaexp 

are allowed to relax for sufficient time, the degree of relaxation to the Maxw(‘llian 
at the corresponding temperature is studied using the nouns defiiu'd in 
eqs. (3.38)-(3.41). These are then referred to as relaxation measun's or norms. 
Such experiments can be used to study numerically the spatial and t<*inpt)- 
ral convergence of various finite elements in the context of the Fokkc'r-Planck 
equations. 

3.4.2 Spatial and Temporal Convergence 

In order to ascertain the performance of various finite elements with the tirin' stejs- 
ping schemes, the results of relaxation of a Gaussian to Maxwellian is pr<'S('nt('d for 
a single species in Table 3.2. All the four measures are shown at lOr (r-relaxation 
time) for various finite elements and time integration schenriis. The rt'sults are 
shown, in percent, for two different total number of grid points (N~()l and 121). 

It is seen that by the time of lOr, the species has almost reached tins (‘(juilibrimn 
state. Since the Crank-Nicolson scheme is second order accurate, its results can bo 
considered to be more accurate (Numerical oscillations were absent in th<^ case of 
a single species). The table shows that backward scheme slows down the <>(iuilibra- 
tion or relaxation process compared to the Crank-Nicolson scheme exc('pt with the 
linear elements with 61 points, the results of which can be considt'red to b<^ th<i U'ast 
accurate. Galerkin scheme is somewhere in between the other two scherncis. Siiute 
it is only conditionally stable, this scheme has not been used in later caU'ulations. 
Among the finite elements, linear elements with 61 grid points show a slightly dif- 
ferent trend. Such trends disappear when the At is changed, as shown lat(^r in 
Fig. 3.3. In general, the linear elements do not show consistent brdiavior with less 
number of grid points. This is explained in the following paragraphs. 
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Table 3.2: Relaxation measures in percent of a Gaussian at lOr with different 
finite elements and time stepping schemes (At = 1/16, L-Linear, Q-Quadratic and 
C- Cubic) 


El. 

S 


N = 

= 61 



iV = 

121 




Loo 

Li 

L 2 

Lrms 

Lqo 

Li 

L 2 

Lrms 

L 

CN 

3.150 

1.053 

1.102 

1.963 

0.448 

0.352 

0.369 

0.413 


GA 

2.713 

0.796 

0.832 

0.160 

0.803 

0.584 

0.618 

0.743 


FI 

1.890 

0.315 

0.334 

0.936 

1.495 

1.030 

1.089 

1.374 

Q 

CN 

1.967 

0.798 

0.847 

1.249 

1.488 

0.796 

0.847 

1.249 


GA 

2.337 

1.029 

1.081 

1.572 

1.859 

1.027 

1.087 

1.481 


FI 

3.033 

1.467 

1.546 

2.183 

2.561 

1.464 

1.544 

2.094 

C 

CN 

1.384 

0.799 

0.850 

1.165 

1.324 

0.796 

0.844 

1.144 


GA 

1.757 

1.030 

1.092 

1.492 

1.697 

1.027 

1.087 

1.470 


FI 

2.459 

1.468 

1.549 

2.108 

2.400 

1.464 

1.544 

2.083 


Convergence with Mesh Refinement 


The spatial behavior of various elements is assessed by studying the relaxation of a 
Gaussian distribution for various total number of grid points keeping the time step 
(At = 1/16) constant. The results are shown in the form of Too and L 2 measures 
at 6t in Fig. 3.2 for different elements. Crank-Nicolson scheme has been used to 
obtain the results shown. It is observed that both quadratic and cubic elements 
closely follow each other. However the linear elements with coarse grid spacing (less 
n nmh pr of grid points) are slow in relaxation compared to the other elements which 
are more accurate. Only on a sufficiently refined grid the relaxation observed using 
the linear elements is same as that obtained by the higher order elements, indicating 
convergence. Similar behavior is seen with the backward scheme also. 


central library 
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Figure 3.2: and L 2 relaxation norms for a Gaussian at fir for tUfFert'ut grid 

spacings 


Convergence with Time Increment 

The temporal behavior of various finite elements is studied in the Kam<‘ manner. 
For a fixed total number of grid points N, time step was varied from 1/2 to 1/64 
T for the Crank-Nicolson scheme. The results are shown in Fig. 3.3 for a singl(^ 
species. Here again the linear elements exhibit a slightly nommiform b<*havior. Both 
the quadratic and cubic elements show good convergence. In general larg<‘r time 
steps slow down the relaxation. As the step size is decreased furtlu'r th(‘ .slop<'s of 
the measures change slowly indicating a steady state like situation. This kiiui of 
behavior is seen when the number of grid points is also vari(>d. It i.s also true with 
the backward implicit scheme. These observations suggest that (piadratic and cubic 
elements are better than the linear elements in terms of couv('rg<‘nce ami ac<'uracy. 

The errors that occur in the case of multispecies calculations have also b('en 
studied. Two Maxwellian species (D-deuterons and T-tritons) with <Hiual (uiergy 
were allowed to evolve with a fixed time step At = 1/16 and the <‘rrors introduced 
purely due to numerical evolution were studied. The velocity cut-offs w(;r<‘ set 
first to a common value for both species and later to 6vtH (multi-velocity) for each 
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Tine step 

Figure 3.3: Loo L 2 relaxation norms for a Gaussian at 6t for different time steps 


species. The errors are shown in Table 3.3 for the second species T. Error norms 
indicate that the multi-velocity scale approach has reduced the errors of the species 
T considerably. This is true with other error norms also. The error norms of D also 
have shown similar improvements, even though in this case the difference between 
the cut-off velocities is not large. The reduction in errors with the higher order 
elements as well as with increasing number of grid points is clearly brought out in 
this table. 

Further to check and validate the implementation of the multi-velocity scale 
model, equipaxtition of energy between deuterons and tritons in a D-T plasma has 
been studied by using both the co mmo n and multi-velocity dimensionless scales. 
Species with different energies at i = 0 {Ed = 75 keV and Er = 15 keV) were 
followed for about 0.4 s as shown in Fig. 3.4 using 30 cubic elements {N = 91). The 
equilibration using a common velocity scale is shown in continuous solid and dotted 
lines. The circles are that of multi-velocity scale model. The results match with 
each other very well. This co nfirms the validity and the proper implementation of 
the multi-velocity scale model. The cumulative error in energy is 0.98% of the total 
energy for the Crank-Nicolson scheme and is 1.9% for the backward scheme. 
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Time (secs) 

Figure 3.4: Equipartition of energy between deuterons and tritons using different 
velocity scales 

3.5 Computational Aspects 

Certain aspects related to density and energy conservation, controlling the time step, 
and changing the velocity scale of a species during the computation are discussed 
in this section. These features have been incorporated in the program to avoid 
cumulative errors in densities and energies and to improve the accuracy of the 
numerical solution as well as to reduce the number of time steps required. 

3.5.1 Density and Energy Conservation 

When the Fokker-Planck equations are solved numerically, it is not necessary that 
both density and energy can be conserved exactly even when large number of grid 
points and small time steps are used. By choosing proper boundary conditions only 
one of them can be conserved (Whitney, 1970). When electrons are modeled as a 
general species, very small time steps are needed. Therefore even small deviations 
in solutions can propagate further as the species evolve to reach steady state. In 
order to maintain conservation either special schemes must be devised, or it may 
be imposed externally. In the present work an attempt has been made to impose 
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the conservation conditions from outside after every time step, or a fixe<i number 
of time steps. 

In order to conserve both density and the total energy ('xactly, it is uecvssary 
to slightly alter the computed distributions using a suitable weighting function for 
each species, i, say + bwi{x) when a, « 1-0, b « 0.0, the comp..t<>d solution is 
expected to approximately conserve the density and energy. Tlu^ piii<mM't,<’rs, a, and 
b, equals the number of constraints to be satisfied. In tins work n'.(x) an' chosen 
to be Gaussian centered around the instantaneous peak of the th<^ spc'c i<'s i. The 
density a n d energy constraints in one dimensions are specifie<l as follows 

f {ai + bwi{x)) fi{x)A'KX^ dx = ni, (3.42) 

Jo 

N ^ 

'^miv'f{ai + bwi{x)) fi{x)2irx*dx = (3.43) 

i=i >-i 

where N is the n umb er of species. The right hand sides of eqs, (3.42) and (3.43) are 
the densities and the total energy to which the computed distribution functions are 
constrained. When source and loss terms are present the cons(^rvation of (k'nsity 
and energy are maintained in the same way by taking into account th<^ rt'spc'ctive 
loss and gain of individual species. 


3.5.2 Controlling the Time Step 

In the simulation of multispecies fusion plasmas it is always desirabh^ to hav<! con- 
trol over the time steps. Adaptive time stepping schemes can be used. In that case 
species are integrated in time with different time steps. SuV>S('(iuently the coeffi- 
cients and the distribution functions have to be interpolated mnong tlu’: specic^s. If 
implemented successfully adaptive schemes are desirable to study even a simi)le re- 
laxation of arbitrary distribution to Maxwellian, where the first few relaxation times 
are more important because by then the spectrum takes the shape of ai>proxiiuate 
Maxwellian. Later on the distributions relax slowly to reach the mutual Maxw<*llian. 
In the present work time steps are controlled in the following ways; 



3.5. COMPUTATIONAL ASPECTS 


57 


i. Constant dimensionless time steps Ar are chosen based on the relaxation time 
of each species which is adjusted if a significant change in the energy of the 
species occurs during the evolution. 


ii. The time steps Ar are controlled by the energy norm of the error in the 
computed solution. 

The number of time steps required in the second option depends on the energy 
accuracy desired. This is elaborated in Chapter 4, where adaptive node distribution 
in (velocity) space is also studied. 


3.5.3 Velocity Scale Refinement 

Multi- velocity scale models together with density and energy conservation, and with 
a strategy for varying the time step, improve the numerical solutions in terms of 
accuracy and speed. The cut-offs of individual species are chosen based on the 
energies at f = 0 and the expected values in the steady state. If reactions and 
other auxiliary heating are present, cut-off velocities are chosen to give adequate 
representation of particles at high energies (tail) as well. 

Even with multi- velocity scale approach, as pointed out in Section 2.3, choosing 
higher cut-offs in the beginning itself to represent species with low energies can in- 
troduce erroneous oscillations making the distribution functions negative, especially 
at the tail end. As the species gain more and more energy due to collisions and 
other mechanisms such negative values can possibly improve slowly. Therefore if 
the velocity scales/cut-offs can be refined during evolution such spurious oscillations 
can be controlled or eliminated. This is particularly important when the electrons 
are treated as a general species. 

A strategy has been used in the present work to refine the velocity scales during 
evolution. If the change in energy of the species exceeds a specified value, given as 
input, then the corresponding velocity normalization constant is modified to a 
new value u"®". Correspondingly the distribution function is projected on to the 
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new grid with the transformations 


(3.44) 


fnew ^ fM (3.45) 

J ^|jnetuJ3 •' ’ 

where f'^‘^ and /”®"' are the distribution functions on the old and u<>w grid.s. Th(^ 
projection on the new grid makes use of interpolation. The only probhun ass()eiHt<Hl 
with such refinement scheme is the projection error. This can be easily controlU'd by 
using higher order interpolation schemes, if at all needed. It is shown in Chapt<'r 4 
that such projection errors are quite negligible. 

Refinement of velocity scale can be better illastrat(^d with an exainpk'. Again 
in a D-T plasma, ions with initial energies Ed = 110 keV and Et -- 10 keV w('r(^ 
allowed to equipart the energy over a period of time irsing; (i) a common wdocity 
cut-off, (ii) different velocity cut-offs resulting in different dimensionU'ss velocity 
scales, and (iii) different cut-offs with velocity scale refimmwuit (eiwagy incrtuiK-nt 
5 keV). The non-negative nature of the solutions at every time st('p is stu<li<‘d by 
summing up the number of grid points at which the distribution function is negativ(5. 
The sum of the negative points is obtained for 0.5 s using 250 st<q)s. Num<Tical 
experiments were carried out for these three cases for various finit(‘ eh'inents and 
for different total number of grid points N. 

The total number of negatives in 250 steps, and the fraction of grid points on 
which the distribution function becomes negative, are shown in Tabl(‘ 3.4 for all the 
three cases mentioned above. The advantages of the multi- v(?lo<'ity scab' model as 
well as the velocity scale refinement are self-evident from this table. In or<ier to 
improve the solutions, in general more number of grid points ar(^ n<‘ed(><i (Chang 
and Cooper, 1970; Harrison, 1988). Such a trend is observed Imre also. 

The energy of the deuterons at 0.5 s is shown in Table 3.5 for various cases. 
It can be seen that the quadratic and cubic elements show acceptable results from 
the view point of energy exchange even with less number of grid points. Furtlu’r, 
it also shows that the velocity scale refinement does not introduce any <*rror in the 
calculated energy of the species. Thus it may be inferred that the velocity tiormal- 
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Table 3.4: Negative values in the spectrum for common and multi-velocity scale 
models (L-Linear, Q-Quadratic, C-Cubic, and N-Number of nodes) 




Common 

Multi-velocity 

With refinement 

El. 

N 

Total 

Fraction 

Total 

Fraction 

Total 

Fraction 

B 

31 

2286 

0.2950 

718 

0.0926 

109 

0.0141 

B 

61 

3575 

0.2344 

349 

0.0229 

21 

0.0014 

B 

91 

4449 

0.1956 

16 

0.0007 

9 

0.0004 

B 

121 

4960 

0.1640 

0 

0.0 

0 

0.0 


11 

1500 

0.1935 

85 

0.0110 

12 

0.0015 


m 

1247 

0.0818 

84 

0.0055 

10 

0.0007 


91 

1387 

0.0610 

0 

0.0 

4 

0.0002 


121 

2127 

0.0703 

0 

0.0 

0 

0.0 


11 

2296 

0.2962 

721 

0.0930 

18 

0.0023 



3654 

0.2396 

529 

0.0347 

10 

0.0007 



4671 

0.2053 

332 

0.0146 

5 

0.0002 



5383 

0.1779 

95 

0.0031 

6 

0.0002 
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Table 3.5: Energy of deuterons at 0.5 s in a D-T plasma obtained by using vjirions 
finite elements 


El. 

N 

Common 

Multi- velocity 

With r<>fin<'m<'nt 

L 

31 

66.78 

62.18 

62.11 


61 

62.30 

61.49 

61.52 


91 

61.76 

61.42 

61.42 


121 

61.59 

61.40 

61.39 

Q 

31 

61.21 

61.39 

61.33 


61 

61.38 

61.37 

61.35 


91 

61.39 

61.37 

61.34 


121 

61.39 

61.37 

61.34 

C 

31 

61.31 

61.34 

61.31 


61 

61.38 

61.37 

61.34 


91 

61.39 

61.37 

61.34 


121 

61.39 

61.37 

61.34 
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ization constants (cut-offs) do contribute to the negative values in the distribution 
functions, and the choice of velocity cut-offs as well eis their refinement during the 
evolution, can play an important role in the accuracy of the computed solutions. 

3.5.4 Structure of the Program 

The programs that have been written here are primarily used to solve the Fokker- 
Planck collision operators based on the finite element method. In addition to the 
features mentioned above, other aspects Hke iterating coefficients over the solutions, 
including source and loss terms, etc. have also been incorporated. The details about 
source and loss terms are given in Chapter 6. The codes can be used to model an 
arbitrary number of general £ind Maxwellian background species, including electrons 
as a general species. The electron population is modified accordingly to satisfy the 
macroscopic neutrality. A flow of the algorithm may be visualized broadly in the 
following sequence: 

1. Read in the input data and other control parameters. 

2. Calculate density, energy, initial distribution functions, etc. 

3. Compute Fokker-Planck coefficients using trapezoidal rule or Hermite scheme 
as specified. 

4. Discretize the equations including source and loss terms and solve the system. 

5. Iterate the coefficients, if desired. 

6. The following options are included: 

• impose density and energy conservation conditions. 

• use velocity scale refinement. 

• change the size of the time step. 

7. Calculate density, energy, and other quantities of interest. 

8. Update the data of source and loss terms, if such terms are included. 
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9. Repeat steps 3 through 8 till the specified exit conditions ar(‘ .satisfi<'d. 


3.6 Finite Difference Solutions 


The purpose of implementing finite difference method of solution for Fokker-Plauck 
equations is twofold. First, the approach of using different dinu'nsionless vtdocity 
scales for different species is tried in this context and compart'd with tht' finite 
element method presented earlier. Secondly, analytical limit of collision opt'rator at 
the boundary derived in Chapter 2 is also implernt'nted and t<'st('<l. 

In the finite difference method, the derivatives in the governing t'tiuations are 
approximated using finite differences over a numerical grid in th<' computational 
domain. Usually the derivatives are approximated with tlu‘ (a>ntral differenct's. 
This approximation translates the differential equation into a system of algebraic 
equations. The 1-D Fokker-Planck equations solved, without any sourc<' and loss 
terms, are given in eq. (2.83). The spatial derivatives are <iiscretiz<'d a<'curate to 
second order. The derivatives are discretized at a node i as 


d(Af) 

dx 


dx [ dx, 


i-^+i fi+i fi-i)/2Axi, 

7i«-/A 


1 

Axi 




AXi^.i/2 


^®»-l/2 


where 


Axi _ -{xi+i ~ Xi_i), 

^®i±l/2 = ± (Xi±i - Xi), 

1 


Bi 


± 1/2 




The above approximations are vaUd for all the inte 


(3.46) 


(3.47) 


(3.48) 

(3.49) 

(3.50) 


interior nodes, except the boundary 
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nodes i — 1 and i = J at x = 0 and x = 1.0 respectively. 


3.6.1 One Step 0-Method of Solution 

The equation (2.83) is advanced in tinae from time step n to n + 1 using the ©- 
method as 


^ = 0£(r^^) + (1 - Q)C{r), (3.51) 

where C{.) is an operator containing the right hand side of eq. (2.83), i.e. 

,df] 




Af + B 


dx 


(3.52) 


Finite difference discretization results in the following tridiagonal system 
tti/i+l + A/» + 

The coefficients and Si are given by 

Oi 

A ^ 

Ti ^ 

6i ■ 

where 


(3.53) 


Oi 


hi 


1 

CD 

(3.54) 

- © hi. 

(3.55) 

cy 

CD 

1 

(3.56) 

(1 - e)(<.i/r« + iiK + 

(3.57) 

1 ] 
x^Aali Y 2 AXi4.i/2/’ 

(3.58) 

1 ( 1 BtiA ] , 1 

x?Axi \^Axi+i/2 Axi_i/2y At’ 

(3.59) 
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l__ (_^ + Billl-] (3.60) 

^ xj^Xi V 2 Axi_i/2 j ’ 

for 2 < i < / - 1. The coefficients A and B are evaluat<'(l using tlu' values of / 
at step n. The above system is solved by using a standard tridiagonal al.-vrithm 
(Richtmyer and Morton, 1967). 

3.6.2 Approximation on the Boundary 

The boundary nodes need special attention, particularly at x — 0 du(* to t lu' sin- 
gularity. The Fokker-Planck collision operator at a; = 0 can h(> treat('d (‘ither 
numerically or analytically. Both the treatments arc iniplenuuited Iwrc. 


Numerical Approximation 

McCoy, Mirin and Killeen (1981) approximated the collision op<'rator at no<i<* i ~ 1 
numerically and the approximation is given by 


1?^ 

x^ dx 


2=0 


X 


^ 3/2 ’ 


3/2 


(3.61) 


where G3/2 is evaluated at the mid point X2/2. The coefficients of th<‘ difh'n'iice 
equation at the boundary can be written in terms of 


ai 


Pi 


-0 


-0 


3 + ^ A^ + Af 

®3/2 V2Ax2_i/2 4 ^ 

3 {B^ + B^ 


''3/2 


V2AX; 


' 2 - 1/2 




At’ 


(3.62) 

(3.63) 


and 


7i =0. 


(3.64) 


The coefficient 61 is 
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Si = (1- e)(aif^ + 4 - (3 65 ) 

The errors introduced due to numerical evolution for a single species are shown in 
Fig. 3.5 along with the finite element results. It can be seen that finite difference 
method is slightly more accurate than the linear elements, whereas the quadratic 
and cubic elements are much better than the finite difference solutions. This clearly 
shows that finite element solutions exhibit better spatial resolution. The behavior is 
the same with all the four error norms. Both Crank-Nicolson and backward schemes 
have shown negligible difference. Further, the figure shows that as the total number 
of nodes are increased from 31 to 121 the improvement in the accuracy of solution in 
the case of finite difference method is slower compared to the finite element method 
whose accuracy improves by almost two orders of magnitude. 


Analytical Approximation 

The equation to be solved at the boundary i = 1 is eq. (2.113), where d^f/Ox^ is 
approximated as 


_ 2(/2 - / i ) 

dx^ (Ax2-i/2)2’ 


(3.66) 


with the condition dffdx = 0. The coefficients ai,/3i, and 71 are given by 


Oi 


-e 


2(7* 


AX3/2’ 


ft = -e ' 


AX: 


3 / 2 , 


At’ 


7i 


0 , 


(3.67) 

(3.68) 

(3.69) 


and 61 is the same as in eq. (3.65). 

The results with the above approximation are shown in Fig. 3.6. The approxi- 
mation based on the analytical limit improves the accuracy (coo) of the numerical 
solution by almost a factor of 2. This trend is seen with Crms also. For a Gaussian 
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30.0 50.0 70 0 90.0 110.0 130.0 30.0 50.0 70.0 90 0 110,0 130.0 

Grid points (N) Grid points (N! 

Figure 3.5: Comparison of numerical errors in. the finite element and finit<> differtuice 


methods 
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30.0 50.0 70.0 90.0 110.0 130.0 30.0 50.0 70.0 90.0 110.0 130.0 

Grid points (N) Grid points (N) 


Figure 3.6: eoo errors of numerical and analytical approximations at the 

boundary for a Maxwellian. 

input, the steady state values are given at 6r in Table 3.6. In this case the analytical 
approximation has increased Too ^2 measures of relaxation. This is correct be- 
cause the numerical approximation at a: = 0 introduces faster relaxation compared 
to the one based on analytical limit of the operator at the boimdary, whose values 
are closer to the finite element relaxation norms, which are shown in Table 3.7 for 
identical input. It can be seen that the norms based on the analytical limi t at the 
boundary are some where between those obtained by the use of quadratic and cubic 
finite elements. 

3.7 Treating Electrons as a General Species 

Errors in the numerical solutions for a single species have shown encouraging perfor- 
mance with the multi-velocity scale approach. Nevertheless, when the ion-electron 
interactions are modeled, the errors caused by the interpolation and integration 
may be significant, due to large difference in the velocities involved. 

To estimate the numerical errors, experiments with two species -ions and elec- 
trons - were conducted for about 6 ion t’s (r-ion relaxation time) with the im- 
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Table 3.6: Comparisoa of relaxation norms at 6r for a Ganssian input u.siug dlihvmt 
approximations at the boundary (CN-Cranh-NifoLsou. FI-Fnlly ImpU.it) 


N 

Numerical (CN) 

Analytical (CN) 

Numeric 

al (FI) 

Analyt 

ical (FI) 


Loo 

Lrms 

ioo 

Lrms 

LiX) 

Lftns 

r 

r 

*^rmg 

31 

1.720 

1.298 

3.385 

3.244 

2.168 

1.852 

4.271 

4.135 

41 

1.427 

1.200 

2.345 

2.225 

2.384 

2.080 

3.350 

3.170 

51 

1.502 

1.308 

2.025 

1.891 

2.522 

2.264 

3.043 

2.857 

61 

1.576 

1.406 

1.899 

1.763 

2.626 

2.382 

2.942 

2.732 

71 

1.639 

1.476 

1.851 

1.710 

2.708 

2.459 

2.916 

2.690 

81 

1.688 

1.526 

1.836 

1.687 

2.772 

2.512 

2.931 

2.(569 

91 

1.729 

1.568 

1.836 

1.678 

2.822 

2.549 

2.975 

2.661 

121 

1.812 

1.626 

1.862 

1.676 

2.924 

2.611 

2.975 

2.660 


Table 3.7: Relaxation norms at 6r for a Gau.ss'ian input asing finite <'l('in<*nts 


Scheme 

N 

Linear 

Quadratic 

Ctibic 



Leo 

Lrms 

A. 

r 

L.u 

T 

Crank- 

31 

9.767 

6.438 

4.278 

2.455 

2.901 

2.057 

Nicolson 

61 

1.105 

0.4411 

2.640 

1.817 

2.063 

1.739 


91 

0.986 

0.8970 

2.295 

1.746 

2.015 

1.721 


121 

1.379 

1.250 

2.171 

1.727 

2.009 

1.716 

Fully 

31 

8.316 

5.115 

5.368 

T 39 O : 

’4l)l() 

3.037 

implicit 

61 

1.086 

0.972 

3.761 

2.801 

3.195 

2.728 


91 

2.087 

1.910 

3.420 

2.731 

3.147 

2.708 


121 

2.505 

2.253 

3.297 

2.712 

3.139 

2.702 
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Figure 3.7: e^o and 62 errors of ions and electrons for a Maxwellian input 

position of density and energy conservation conditions. The results are shown in 
Fig. 3.7 for quadratic elements. Error norms Coo and 62 indicate that the errors are 
under control and that the multi-velocity scale model is equally well applicable to 
model the electrons 8 is a general species. 

It may be noted in Fig. 3.7 that the magnitude of errors of the electrons which 
have larger cut-off is less than that of ions. This is due to the representation of 
species on different velocity grids. When the Fokker-Planck coefficients are calcu- 
lated for electrons, the contributions from the ions to each of the points on the 
electron grid can be obtained easily. On the other hand the converse is not true 
because the contributions from electrons to almost all the points on the ion grid 
essentially come from the first few points of the electron grid. This is due to the 
factor of 40-60 in the velocity cut-offs of electrons and ions. This might give rise 
to some error in the ion equations. However as the total number of grid points 
is increased, the difference between the errors of ions and electrons become less, 
especially in the L 2 error. 

Further to study the energy transfer between ions and electrons, a deuterium- 
electron plasma has been modeled. The Fokker-Planck equations were solved both 
by the finite element and the finite difference methods. When ions and electrons me 
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at low temperatures, in the absence of external influences, .be e.n.r«y transfer (r. 
luxation) can be predict by the Spitzer model approxuna.ely (bp.tr.er r., 1%2). 
At low temperatures it may be assumed tlmt the dlstnbutinus renuun Maxwelhan 
during relaxation. However even if one of the species is at high temperature, the 
distributions will not remain Maxwellian during relaxation. This is beca.tse the 
particles in the low energy region gain more energy front the h.,t .species causing a 
distortion in the distribution function (Killeen, Heckmtte att.i Boer, l-)l,2). 

Energy transfer between ions and electrons at low temperattne has bee.t sttnlusl 
mainly with a view to compare the Spitzer model with the Fokker-Ph.nck results, 
using the methods developed in this work. Acconling to Spitzer, in a Maswellkui 
plasma the equipartition of temperature from ions to electrons is governed by 


dPe I Td- Te 

Hi "'to {T,+ {mJmn)TnW^ 


where the parameter to is given by 


= 1.8 X 10' 


mo 


(3.70) 


(3.71) 


and the symbols have their usual meanings. The temp('ratur('.s are in unit.s of 
keV and the density in Equation (3.70) can be int<'grat(‘d <«iv.sily by using 
Te + Tb = To = (constant). The solution is (Prudnikov, Brychkov and Mari<-hev, 
1986) 


t = to 


(0) - g^l^) + e {g''\0) - g'l^) + ^ in 1 ! 


(3.72) 


with a = mJmD, b = —Tol2, c = {1 — a), d = aTo, e = {d - In:), x - f- b, and 
g = {cx + e). $(f) is given by 


^ (cx(t) + e)^/^ - 
(cx(t)+e)V2+ei/2’ 


(3.73) 


where $(0) is the value of $ at t = 0. The energy transfer from ioas to <'lectrons 
has been studied for a range of T^/Te = 10 to T^/Tf = 100. For the case of 
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Td ITe = 100, both the finite element and the finite difference results are shown in 
Fig. 3.8. The initial temperatures are Td = 4.0 and T^ = 0.04 keV. 

The results have been obtained using only 61 grid points. Quadratic and cubic 
elements have produced the same results. The multi-velocity scale model hais made 
this possible which otherwise would have needed at least 150 points (Whitney, 
1970). The figure shows that the Spitzer model is in good agreement with the 
Fokker-Planck calculations. The maximum difference between the two models is 
about 6%. In this case both the finite element and the finite difference solutions 
agree with each other very well. 

Energy transfer from electrons to deuterons also has been studied in the same 
manner. Here the Spitzer model is given by 


dTp^^l T,-Td 

dt ^ to (Te + {me/mD)TDyf‘^' 


(3.74) 


The parameter to is given in eq. (3.71). The solution is 


t = to 


(0) - p"/') + b - g^f^) + 


^ $( 0 ) 

2 '$(t)'J 


(3.75) 


with a = me/mD, 5 = (1 + a)To/2, c = (1 — a)/2, x — To — 2Tb, and g — {cx -f h). 
$(t) is given by 


_ (cx(t) + 

(ca:(t) -t- 5)^/2 + {, 1 / 2 - 


(3.76) 


Here also the temperature ratio Tp/Tg was varied from 0.1 to 0.01. Keeping other 
inputs the same, the result of a test case is shown in Fig. 3.9 for the initial tem- 
peratures Tp = 0.04 keV and T^ = 4.0 keV. In this case too the Spitzer model is in 
very good agreement with the results of the finite element method. The maximum 
difference in the finite difference case is about 8%. Therefore it can be concluded 
that the Spitzer model predicts the energy relaxation under Maxwelhan conditions 
adequately. 

Further a study of energy transfer from hot ions to cold electrons and vice versa 
has also been made using the multi-velocity scale Fokker-Planck equations. Ions 
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Time (secs) 

a) Finite element method 



Figure 3.9: Energy transfer from electrons to ions at low temperatures 
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a) Deuterons to electrons b) Electroiw to 

Figure 3.10; Energy transfer between ions and electrons Finite (‘lement method 


with 105 keV and electrons with an initial energy of 15 w('r<' alIow<‘d to <'<nupart 

the energy for about 5 s. Total energy and density ar«> conserved in all tin' cases 
presented here. The results are presented in Figs. 3.10 and 3.11 for the finite (‘huiK'nt 
and the finite difference methods, respectively. In gimeral both th(' nwt Inals closely 
follow each other. However, in the case of the finite diff(‘ren(‘<> method, to obtain 
the results similar to the finite element method, the total numln'r of grid points had 
to be increased from 61 to 121. 

The studies have shown that to improve the accuracy of the finiti* <litference 
solutions in general more number of points are n<‘<'ded. On the other hand the finite 
element method is able to produce accurate results with h\ss numb<>r of I'h'un'uts. 
The CPU time of the various finite elements and the finite differeuci* methoii are 
given in Table 3.8. The run times given are on the CONVEX-22() systinu with 01 
level optimization. CPU times are shown for optioas (I) and (11). In (II), th<‘ eimgy 
and density constraints as explained in Section 3.5.1 are imposcai. It shows that the 
imposition of these conditions, which requires the evaluation of a few int«*grals at 
each time step, consumes a reasonable fraction of the total tim<>. 
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Time (secs) Time (secs) 

a) Deuterons to electrons b) Electrons to deuterons 

Figure 3.11: Energy transfer between ions and electrons-Finite diJfference method 

It can be seen that although the finite element method takes about twice the 
time as compared to the finite difference method for the same number of nodes, the 
improvement in accuracy is much more. This can also be seen from Fig. 3.5. When 
the number of nodes are increased from 31 to 121 the errors in the finite difference 
method are reduced by less than a factor of 10, while the finite element method 
shows an improvement by almost two orders of magnitude. This indicates that the 
spatial resolution of the numerical solution of the Fokker-Planck equations based 
on the finite element method is better than the finite difference method. 

3.8 Summary 

The experience in this chapter has shown that the solution of multi-velocity scale 
Fokker-Planck equations based on the finite element method can be very useful in 
the area of numerical modeling of fusion plasmas. It is found that for the same com- 
putational effort the finite element method performs better than the finite difference 
method. The work reported in this chapter can be summarized as follows; 
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Table 3.8: CPU times in minutes for the hnite element (FE) and the Unite ,lifr,T..nre 
(FD) methods 


Option 

N 

Linear 

Quadratic 

Cubic 

1 FI) 

I 

31 

2.3 

2.1 

2.3 

l.O 


61 

4.9 

4.4 

4.9 

2.3 


91 

7.8 

7.0 

7.8 

3.9 


121 

11.0 

10.0 

11.0 

U,. . , . 

5.9 

II 

31 

3.0 

2.8 

3.1 

i 

1.7 


61 

6.3 

5.8 

6.3 

3.7 


91 

9.9 

9.1 

9.9 

6.0 


121 

13.8 

12.8 

13.8 

8.6 


i. Multi-velocity scale Fokker-Planck equations in 1-D have Ixh'u s<>1v<h1 by both 
the finite element and finite difference methods. 

ii. Strategies have been presented to vary the time steps and also t.o sat isfy the 
conservation of both density and energy at cv(’!ry st(^p. 

iii. It is found that the choices of velocity nonnalization constants (cut-otfs) do 
improve the non-negative nature of the numerical solutions significantly. Fur- 
ther, varying the cut-offs during the course of evolution of plasma also reduces 
the problem of solutions becoming negative at tlu' tail. 

iv. Quadratic and cubic finite elements are better than the linear <'lem<'nts in 
terms of convergence and computational accuracy. 

V. Electrons have been treated as a general speci(« togeth<*r with ions. Using 
the multi-velocity scale model, the energy traasfer betwenm ions and eh'ctrons 
has been studied using same number of grid point.s. Th<' results of Spitzer 
model which is valid under Maxwellian approximation havt' been found to be 
in good agreement with the Fokker-Planck solutions. 



Chapter 4 


Discretization Errors and Optimal 

Distribution of Nodes 


4.1 Introduction 

When a dimensionless velocity scale characteristic to each species is used in the 
simulation of multispecies plasmas, the relevant velocity ranges are better modeled 
even with a uni form distribution of nodes. With this approach, the set of nodes for 
the species with the largest velocity cut-off, say electrons, need not be the union 
of nodes for other species plus additional nodes to cover the remaining part of the 
electron velocity. The number of nodes for each species can be chosen independently, 
and if desired, the same number of nodes can be used for all the species. The 
nodes, however, do not correspond to the same physical velocities, and appropriate 
interpolations are required for calculating the Fokker-Planck coefficients. The fact 
that this approach can be successfully exploited both in the finite element method 
and the finite difference method has been demonstrated in the previous two chapters. 

Distribution functions, however, vary by a few orders of magnitude even when 
a dimensionless velocity scale most appropriate for each species is considered. Con- 
sequently, the discretization errors associated with a uniform distribution of nodes 
show a large variation. In this chapter, the local discretization/truncation errors 
are studied for typical distribution functions encountered in the simulation of fusion 
plasmas. These include the equilibrium Maxwellian distributions and the Gaussian 
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distributions which can be used for the characterization of fast ions or reaction 
products. Redistribution of nodes with a view to immmi/.e tlie errors and aciiu've 
a more uniform error distribution is discussed. When* po.s.sil,le, exact errors iue 
. compared with the computed errors. Errors for ions and (‘lectrous are compared, 
both on a uniform grid as well as on an optimal grid. An adaptiv,* time stepping 
scheme based on the norm of the energy error, as mentioned in Hecti<vn 3.5.2 is also 

implemented. 


4.2 Error Estimates 

The Fokker-Planck equation used in this chapter is in the form of t'cp (2. 104): 


dr 


df d^f 


£(/). 


(. 1 . 1 ) 

(4,2) 


where A*,B*, and C* are as in Chapter 2 (eqs. 2.105 - 2.107), and C denot<*s the. 
operator on the right hand side. Three different methods to estimate the nodal 
errors were considered; 


i. Developing the truncation error term for the operator C aualyti<'ally in t(‘rms 
of the higher order derivatives, 

ii. Use of higher-order finite difference formuUie in conjunction with th<^ c(‘utral 
difference formulae for approximating the operator C and trc'ating th<' <iifFer- 
ence as the error estimate, and 

iii. Approximating the operator £ on a fine and a coarse gri<l, and tr<'ating the 
difference as the error estimate. 

The first method requires the computation of the higher order (h'rivatives. Further- 
more, analytical expressions for the leading terms in the truncation (‘rror(s) become 
intractable on a non-uniform grid except for the standard three-point central dif- 
ference schemes. Comparison of error estimates with exact errors, which can be 
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done at t = 0 (see Section 4.4) has shown that generally the estimates obtained by 
method (iii) are more accurate as compared to those obtained by the methods (i) 
and (ii). Method (iii) is also the simplest to use. It is this method which has been 
used in this work, and can be expleiined as follows. Let Cax and C 2 Ax be the two 
approximations of the operator C on the fine and the coeirse grids respectively. If 
€ax and € 2 Ax are the corresponding discretization/truncation errors, and n the order 
of convergence for the discretization scheme, then 

= EAxif) "b ^Ax) 

= .^2Ax(/) ^lAxi 

= CAa:”, (4-5) 

e2Ax = C'(2A®)" = 2"eAx. (4-6) 


(4.3) 

(4.4) 


dr 


where 




Subtracting eq. (4.3) from eq. (4.4), and using eq. (4.6), 

1 Eax ~ F-iax 1 


I^Ail 


2 "- 1 


(4.7) 


For obtaining E^Ax-, alternate points are skipped from the fine grid. 


4.3 Redistribution of Nodes 

The purpose of redistribution of nodes is the minimization of and a more uni- 
form distribution of the discretization errors. The optimal distribution is obtained 
by iteratively adjusting the grid spacing in inverse proportion to where n is 
the order of convergence. Since the nodes in the present work are distributed on 
the dimensionless velocity scale extending from x = 0 to 1.0, the relation giving the 
new nodal distribution (spacing) Ax^ in terms of the old distribution Ax^ can be 


written as 
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Ax' 


E/.iAxi/|£i|'/«- 


(4.8) 


Starting from a uniform grid at t = 0, generally 4 to 5 iterati<>n.s l)as(‘<l on «*({. (4.K) 
are required to minimize €max- At later time steps oik; it«‘nit ion, after a fixed number 
of time steps or time of integration is found to bo a<ie«}uat<'. If the I'rror f, for a 
particular node i is below a certain cut-off limit, it is this limit which is us«’d to 
avoid excessively large grid spacing around such node's. 

Figures 4.1, 4.2, and 4.3 present the results obtaim'd with 51 and 101 node.s, 
using both the standard as well as higher order five-point difh'ri'm-e formulae basi'd 
on Lagrangian polynomials. Figures 4.1 and 4.2 show the ilistribution of relative 
nodal errors {ci/fr^ax) on the uniform and the optimal grid for a Maxwelli.m distri- 
bution function. It is seen that the maximum discri'tization ('rror on the optimal 
grid is reduced by a factor of about 7 in the case of thr«'i'-point approximation, 
and by more than an order of magnitude in the case of fivi'-point aj>i)roximaf ion. 
The results are similar for a Gaussian distribution wliich an' shown in Fig. 4.3. It 
should be mentioned here that the error distribution shows a number of maxima 
.and minima both on the uniform and the optimal grid. This is due to tlu* highly 
complex and nonlinear nature of eq. (4.1) where the coeffiemnts B\ ami C* in- 
volve various derivatives of the Rosenbluth potentials. Further tlm truncation <'rror 
for the second and third terms in eq. (4.1) combine differently at difh'rent points 
on t e grid, leading to addition or cancellation. More iterations of <«q. (4.K) gem'r- 
ally do not reduce such fluctuations even on the optimal grid. Fig. 4.4 shows the 

ac ua istnbution of nodes after five iterations of eq. (4.H) for a Maxw,-lliaa and a 
Gaussian distribution. 


4.4 Exact and Computed Errors 

The redistributioa of nodes shown in Figs. 4.1-4.3 has b.H.n based on the error 

distributions ex! t J ^ for uniform as well as non-uniform node 

^tnbntrons, exa« ^crefsation errors at t = 0 for known dhstribution function. 
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Dimensionless Velocity Dimensionless Velocity 

a) 51 nodes b) 101 nodes 

Figure 4.1: Comparison of discretization errors for a Maxwellian distribution on 
uniform and optimal grids using standard three-point central difference formula 



a) 51 nodes b) 101 nodes 

Figure 4.2: Comparison of discretization errors for a Maxwellian distribution on 
uniform and optimal grids using a five-point difference formula based on Lagrangian 
polynomials 
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0.0 0 0 

Difnensionleoo VoliuMt / 

a) Three-point central differem^e formula 

-2 



n 0 

Oiireiisiur-ln.V'. '-Ailfir i i.y 


b) Five-point difference formula 

Figure 4 3: Discretkation errors for a Gaussian distribution on a uniform and «l>- 

timal grid with 101 nodes 
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0.0 0.2 04 0.6 0.8 1.0 

Dimensionless velocity 


a) Maxwellian 



b) Gaussian 


Figure 4.4: Optimal distribution of nodes for two different distribution functions 
(101 nodes) 
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-2 



uTi 

fV 

ai 



b) Optimal gri(i 


Figure 4.5: Comparison of exact and estimated errors for a Maxw<*llian distribution 
using central differences 


were computed. For the Maxwellian distribution £(/) = 0, and th<>rof<)ro 1 £ax(/)| 
itself can be treated as the exact error. For a Gaussian distribution, t.he exju:t 
values of C{f) can be obtained using analytical integrations. This how<'v<T is <iuite 
involved. Instead, a 12-point Gauss quadrature together with analytical valiu's of 
the derivatives were used to calculate exact C{f). Comi)arisous with C^j,{f) yields 

the exact ca*. This is then compared with the estimates obtained by using the fine 
and coarse grids. 

The results are shown in Figs. 4.5, 4.6, 4.7, and 4.8. It can be seen that the 
exact errors compare quite well with the estimated errors. In these fig»ir(*s, the two 
errors have been normalized to agree at the point (node) of maximum «‘rror. In the 
distribution of nodes, it is the relative variation of the error along the nod<’s which 
is of consequence. The good agreement between the exact and the (;stimat('d terrors 
inicates that the estimated errors can be used effectively for the purpose of node 

ution to obtain an optimal grid initially, and for changing it at subsecpwuit 
steps during the evolution. 

Figure 4.9 shows the variation of L 


’max error, i.e, Cmax, for diffenmt average grid 




Relative Error Relative Error 
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Diaensionless Velocity Dinensionless Velocity 

a) Uniform grid b) Optimal grid 


Figure 4.6: Exact and estimated errors for a Maxwellian distribution using five-point 
difference formula 



a) Uniform grid b) Optimal grid 

Figure 4.7: Comparison of exact and estimated errors for a Gaussian distribution 
using central differences 
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Figure 4.8: Exact and estimated errors 
order difference formulae 


10 L 



0.0 I) 10 

b) Optimal grid 


for a Gaussian distribution using higlKir 


spacings on the uniform and optimal grids. The reduction in discretizat ion (*rrors 
achieved by an optimal distribution of nodes is clearly brought out in this figur<' for 
all values of the grid spacing. The order of convergcmce n both for th<‘ thr<'<’-point 
and five-point finite difference formulae can be estimated from this figure*. It is 
found to be 2 for the first case and varies between 3 and 4 for th<^ se'cond <'iise. 
It may be mentioned here that the number of iterations reepured bfise'd on <*< 1 . (4.8) 
are relatively insensitive to the exact value of n. Further, the e'xact value* of n is not 
required in eq. (4.7), as only the relative distribution of the errors is of conseepumce. 
Finally, Fig. 4.10 compares the errors for a Maxwellian distribution with a Uurge 
velocity cut-off, which is sometimes needed to represent a bac'kground cold phisma 
species. It is seen that the reduction in €max is significantly large^r (more than two 
orders of magnitude) in this case if an optimal grid is used. 
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-1 



a) For a Maxwellian distribution 


0 



b) For a Gaussian distribution 

Figure 4.9: Variation of error with grid spacing: (U) Uniform grid. (O) Optimal 
grid 
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Figure 4.10: Comparison of discretization errors for a Maxw('lliaa distribution with 
a large velocity cut-off 


4.5 Comparison of Errors for Ions and Electrons 

As has been mentioned in the earlier chapters, and in the introduction to ttu‘ present 
chapter, a characteristic dimensionless velocity scale, x, <'xt«'uding from 0.0 to 1.0, 
has been used in this work to represent all plasma species. Sinc(‘ the physi<'al cut- 
off velocities for ions and electrons are very different, the discretization errors for 
the two species were calculated separately along the lin(\s alrt'ady dis<*u.ssed in the 
previous sections. The same number of nodes were uschI for both th<; spe('i<‘s. The 
calculations were done using both the three-point central differ<'nc(‘s its well as a 
five-point Lagrangian interpolation formula. 

The results are shown in Fig. 4.11. It can be seen that if a dimensionless velocity 
scale is used, the discretization errors are quite comparable for ions and (‘U:ctrons 
using the same number of nodes. This is true both on the uniform grid as well as 
on the optimal grid. The reduction in €max by a factor of 7 to 12 in going from a 
uniform to an optimal grid can be easily observed. 

Figure 4.12 shows changing distribution of nodes as a Gaussian distribution 




RelaU'ie Error 


electrons 
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evolves to a Maxwellian at different times. The number of nodes used is 101, and the 
grid is adapted after every three time steps. The gaps between the dots along the 
x-axis represent the grid spacing at the various stages of relaxation (t is th(^ i(dax~ 
ation time). This figure indicates that the error estimates and the method for grid 
adaptation presented here can be used for studying the plasma evolution starting 
from a highly peaked distribution to a near-Maxwellian equilibrium distribution. 


4.6 Interpolation and Quadrature Errors 

In grid adaptation, as the grid is changed from the old to the new, the distribu- 
tion functions as well as the coefficients in the Fokker-Planck equations have to be 
obtained on the new grid. Since this has to be done a large number of tim(^s, the 
projection errors must be kept small. In order to project the values from one grid 
to other, two points on the old grid closest to the desired point on the new grid are 
located for the purpose of interpolation. Both linear and cubic (Hermitaiu and La- 
grangian) polynomials have been used for this purpose. The errors involved in one 
such step are shown in Table 4.1. It can be seen that the maximiun error per step 
is 0.1% if linear interpolation is used. These errors can be reduced to a negligible 
level if cubic polynomials are used. It should be mentioned here that when different 
velocity scales are used for different species, such interpolations are required at each 
step for the computation of the Fokker-Planck coefficients. However, since these 
coefficients involve the integrals, errors are generally small. These are shown in 
Table 4.2 for 31 grid points. In most of the calculations, the number of grid points 
used is larger than 31 and the actual errors are even smaller than shown in the 
table. 


4,7 Adaptive Time Stepping Based on Energy 
Error 

As mentioned in Section 3.5, a numerical solution of the Fokker-Planck equation 
does not conserve densities and energy exactly. This is due to the discretization 
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Figure 4.12: Optimal distribution of nodes for the evolution of a Gaussian at dif- 
ferent times 
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errors in the velocity space as well as due to a finite time step used in time marching. 
The velocity space errors can be reduced by using a larger number of grid points 
(see Fig. 4.9). In order to control the errors due to time stepping, an adaptive time 
stepping scheme based on the norm of the energy error in the computed solution 
has been used. This is done in addition to the renormalization of the densities 
and energy as explained in Section 3.5. This controls the energy error in each 
step, avoids the cumulative errors in density and energy, and generally reduces the 
number of time steps required. 


Table 4.3 shows the density and energy errors for a single species at each time 
step. It can be seen that the density is almost exactly conserved, while the maximum 
energy error is approximately 0.4%. The corresponding renormalization factors a 
and b based on eqs. (3.42)- (3.43) are also shown in the table. It can be seen that for a 
single species 51 grid points are adequate, while for a multispecies plasma simulation 
involving different types of distribution functions, generally a larger number of grid 
points are required. The error also depends on the time step At. While the energy 
error per step is small, unless renormalized to the desired value, the cumulative error 
in, say, 100 steps can become significant. The table also shows that, eis expected, 
a w 1.0, b = 0.0. 


The energy error shown in Table 4.3 has been used to control the time step size. 
Table 4.4 shows the number of time steps required to evolve a Gaussian species for 
1 s, starting with At = 0.1. If a constant step size is used, 287 steps are required. 
The number of steps required in the adaptive time stepping depends on the energy 
error, AE/E, tolerated per step. For AE/E = 0.5 to 1.0 %, the number of steps 
required is in the range of 40 to 50, while for AE/E — 0.05 to 0.1%, the steps 
required are between 400 to 500. For other values of AE j E, the number of steps 
required are given in the table. It is seen that if reasonable lim its for energy error 
are specified, the number of steps required and therefore computational time can 
be significantly reduced. 
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Table 4.3: Density and energy errors per time step and the nmonuali/ation con- 
stants a and b 


N 

At 

Density error % 

Energy error % 

a 

h 

51 

0.05 

0.00001 

0.0028 

0.9999 

0.00016 


0.10 

0.00002 

0.0115 

0.9996 

0.00070 


0.20 

0.00006 

0.0521 

0.9980 

0.00320 


0.50 

0.00037 

0.3829 

0.9854 

0.023‘K) 

101 

0.05 

0.0 

0.0028 

0.9999 

0.00016 


0.10 

0.0 

0.0115 

0.9996 

0.00070 


0.20 

0.0 

0.0521 

0.9980 

0.00320 


0.50 

0.00033 

0.3774 

0.9856 

0.02390 

151 

0.05 

0.0 

0.0028 

0.9999 

oTooirn" 


0.10 

0.0 

0.0115 

0.9996 

0.00070 


0.20 

0.0 

0.0517 

0.9980 

0.00320 


0.50 

0.00030 

0.3791 

0.9855 

0.023()5 


Table 4.4: Number of time steps required for the evolution of a Gaussian species 
starting with At = 0.1 for = 1.0 s 


Number of 

grid 

points 

r Constant 
At 

Adaptive At 

AE/E 

0.5-1.0% 

0.5-0.25% 

0.25-0.125% 

0.10-0.05% 

31 


52 

103 

189 

512 

51 


50 

102 

183 

506 

71 


46 

89 

200 

470 

101 

287 

44 

93 

163 

445 

151 


44 

81 

161 

416 

201 


41 

79 

168 

403 
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4.8 Summary 

In this chapter the discretization errors in the velocity space have been studied. 
A strategy to redistribute the nodes iteratively to minimize the errors has been 
described and implemented. The maximum error can be reduced by an order of 
magnitude or more if an optimal grid is used instead of a uniform one. The error 
estimates used in the grid adaptation are checked against the exact errors, where 
possible. The convergence of Cmai with the number of grid points has also been 
presented. The discretization errors for ions and electrons have been compared. 
Finally, an adaptive time stepping scheme based on the energy error in the computed 
solution has been presented. 



Chapter 5 


Numerical Solution of 2-D 
Fokker-Planck Equations 


5.1 Introduction 

Numerical simulations using one dimensional analogue of a complex physical sit- 
uation are naturally useful and make the problem tractable. Such models also 
reveal the intricacies involved in the process of simulation. Once the 1-D models 
are developed and implemented successfully, the next logical step in the course of 
the s im ulation studies is to extend the numerical models to the next higher di- 
mension. This is essential to obtain more realistic answers to the actual physical 
situations. Prom the computational point of view, such answers should be obtained 
at an acceptable effort. This gives an opportunity to devise new numerical schemes, 
in addition to improving the existing ones wherever possible, that converge to the 
correct answers more rapidly. 

Numerical simulations of 2-D Fokker-Planck equations are presented in this 
chapter. The equations are solved by both the Galerkin finite element method and 
the finite difference method. Five types of quadrilateral elements namely, linear, 
and the complete and incomplete quadratic and cubic finite elements, are used. 
Discretization of the spatial part of the Fokker-Planck equation is given in detail. 
The equations are advanced in time by the one step 0 method. In the case of the 
finite difference method the equations are integrated by the fully impUcit and the 
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alternating direction implicit methods. As in the previous chapters, a dimension- 
less velocity scale and a dimensionless time most appropriate for each species is 
used. Further, the operators used are as obtained in Chapter 2, that is, without 
singularities at a: = 0 and 0 = 0 and tt. 


5.2 Finite Element Solution in 2-D 

Galerkin finite element formulation of the 2-D Fokker-Planck equatioas in (v, 0) 
geometry proceeds on similar steps as in the case of 1-D fonnulations. Here the 
computational domain is a semi circle with radius Xmax (which is equal to 1.0), aud 
the pitch angle 0 is between 0 and tt. The domain is divided into NgNo curved 
quadrilateral elements, where N^ and Ng are the number of elements in x and 0 
directions respectively. A typical domain and the mesh partitioning is shown in 
Fig. 5.1. The elements that have been tried out in this work are also illastrated 
in Fig. 5.1. The nodes in each element are numbered as shown in the same figure. 
The incomplete quadratic and cubic elements do not have interior nodes. 

It is possible to use either triangular elements completely or a combinat ion of 
both triangular and rectangular elements to discretize the domain. In this work only 
quadnlateral elements are used. The curved elements are mapped into quadrilateral 
elements using isoparametric transformations. 

5.2.1 Formulation 

The 2-D Fokker-Planck equation to be solved for a single species is of the form 

g/ ^ /igG 1 dH\ 
dr dx x^sin.0 d0 J ^ 


where 


(5.1) 
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Finite Element Mesh 



Linear Quadratic (C) Quadrotic (I) 

(4 nodes) (9 nodes) (8 nodes) 



Cubic (C) Cubic (I) 

(16 nodes) (12 

Figure 5.1: Finite element computational domain and various types of quadrilateral 
finite elements 
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The symbols have their usual meanings as given in ChapttT 2. The boundary 
conditions are as given in eqs. (2.69)-(2.72). The residue R is 


df (idG 1 dH\ 
dr dx x^sin0 89 ) 

and its integral is minimized as 


-Qf-s, 


WiR2'irx^sm9dBdx = 0, i=l,. 


where the factor 27rx^sin0 is the 2-D differential area element in radial 
The approximate solution / at a time t is written over an element c as 


(5.4) 


(5.5) 


g(H)metry. 


r 


E Mt) 


i=i 


(5.6) 


where the shape functions are functions of x and 9, and fj, the nodal values 
of the distribution function /, are functions of time r. The number of nodes per 
element is m. Choosing Wi = <f>i and substituting eq. (5.4) in eq. (5.5), 


f r \?l - ( n- 

^ x^sin^ do J ^ 


(f>ix'^ sm9(i9dx = 0, 


i = 1, . . . , m. 


(5.7) 


Integration of the above equation by parts for an element e leads to 
// J sm9d9dx 

"*”/ / ^iQ sin9d9dx 

-J j <f>iSx^ sm9d9dx- f"'U,Osm9f" d9 

^ J X« 

i=i m (5.8) 
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where Xf, Xf+i and 9e, ©c+i are the global x and 6 coordinates of the corner nodes of 
the element e. Substitution of eq. (5.6) in the above equation leads to 



E[(// Sin0d0da:j^ + (7i+/2 + /3 + /4+/5 + l6+/7)/; 


= e[/ f^.Sx^Aueded^+f-"[^fismeY"' + 

■ 


* = 1, . . . ,m 

(5.9) 

where 



/i 


(5.10) 

I 2 

= / f ^A(i>j sine dedx, 

(5.11) 

h 

^ JI ae 

(5.12) 



(5.13) 

h 


(5.14) 

ig 

- //S'S"-* 

(5.15) 

h 

= — J J (f>iQ(f>jX^ sinOdedx. 

(6.16) 


By using the boundary conditions eqs. (2.69)-(2.72) and the values of the coefficients 
A through F at the boundaries, it can be shown that 

t%,H\x = 0, ( 5 - 17 ) 

Jo 0 


because JET = 0 at 0 = 0 and tt. Shnilarly 
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f (biG sin ( 
Jo 


1.0 


d9 = 0, 


(5.18) 


since G is zero at x = 0 and / = 0 at x = 1. The equations (5.9) can tiu'n be 
written as 


t— + uf(T) = v(r), 
dr 

(5.19) 

where 


tij — J j <f)i(l>jX^ sia 9 d9dx, 

(5.20) 

Uij = Ii + 12 Iz F li 15 F Ib F It, 

(5.21) 

Vi = J J <j)iSx^ sm9d9 dx + J G sin 9 

d9 


( 5 . 22 ) 


Element matrices are generated for all the elements in the mesh. These matrices 
are assembled to form the global finite element discretized equations (Dhatt and 
Touzot, 1985) given by 

T| + Uf=V. (5.23) 

Equation (5.23) represents a system of coupled ordinary differential equations. This 
system is solved using the one step 0 -method (time marching scheme). 


5.2.2 Evaluation of Element Matrices and Coefficients 

Each of the elements shown in Fig. 5.1 is represented by isoparametric rectangular 
elements with local variables (^, 17 ) in the x and d directions. The transformation 
of global coordinates (x, 9) to local coordinates (^, t]) can be written in terms of the 
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shape functions as 

m 

i=l 

m 

0 = 


(5.24) 

(5.25) 


1=1 


where Xi and Oi are global coordinates. With the help of above transformations, 
the integrals of element matrices in eqs. (5.20)-(5.22) can be evaluated using Gauss 
quadrature. The differential element area dx dd is written as 


dxd9 = \3\ didr}, 

where |J| is the determinant of the Jacobian matrix 


dx d9 


(5.26) 


dx d9 

di] 


(5.27) 


The elements of J can be obtained by differentiating eqs.(5.24) and (5.25). For 
rectangular elements |Jl is the area of the element, which is known. The derivatives 
d(l>/dx and d(l>/d9 are calculated using the following transformation 


(5.28) 


Using equations (5.24) and (5.28), the integrals in eqs.(5.20)-(5.22) can be rewritten 
in terms of the variables {i,r)) as 


r^i 

dx 

1 

d9 

d9 1 



W J 

"IJI 

dx 

~di^ 



dj) 

. ^ . 


tij = <l>i (t>j 1 J1 sinOdi dr], 

Uij = Jj + iiz + "I" ^4 "b .fs "b -^5 "b 


(5.29) 

(5.30) 
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where 

r+l r+1 Q(h‘ 

^ J-i J-I sine d^dri, 

y+1 f+idSi, 

"" I-i I-I 

^ ~/i /j Sin0<z^£i»7, 

The 2-D integrals in the above equations are evaluated by Gauss 
general 2-D integral is evaluated as 


sia $ 


X,|| 




dr) 


(5.31) 


(5.32) 

(5.33) 

(5.34) 

(5.35) 

(5.36) 

(5.37) 

(5.38) 
quadrature. A 


y+i f+i A', 

i-i 7-1 = 


iV« AT- 


»=1 i=l 


(5.39) 


w ere fe,,,) are Gauss poiuts and »„u,^ are the corresponding weights. It is not 
necessary that same number of Gauss points N, should be used for both directions. 
An easy way to generate two dimensional basis functions is to form the product 

“ ‘-e tor complete elements. However 
r incomplete elements without interior nodes, these fonctions differ slightly. The 

.tape functions Mi, V) for five types ofelements considered in this work, vis., Lear, 
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complete and incomplete quadratic and cubic elements are given in Appendix B. 

The coefficients A through F appearing in the element matrices are approxi- 
mated as described in Section 3.2.2. Their nodal values are interpolated for different 
values of r)j). The other coefficients Q and S are also handled in the same 
manner. 

The elements at a; = 0 require special mention. Here, collapsed rectangular 
elements are used for mapping. In this case, one side of the rectangular elements 
is collapsed onto a point by setting the global coordinates to be the same for all 
the nodes on that side. This approach is followed in order to use the same type of 
elements in the entire domain. 


5.3 Accuracy and Numerical Results 


5.3.1 Error Indicators and Relaxation Norms 

The convergence and accuracy of the finite element solutions can be improved by 
using higher order elements, increasing number of elements and/or by decreasing 
the time step, etc. It is always desirable to obtain accurate numerical result at 
a minimum computational effort. The quality of the numerical solution of the 
2-D Fokker-Planck equation is measured using the error norms as defined in Sec- 
tion 3.4.1. For 2-D these error indicators in percent can be written as 


ei 


62 


|A/. 


Jhllsm. X 100, 

fmax 


f/IA/l x'^sinB d9dx ^ 
/ / sin.9d9dx 


/ / |A/r sinOdBdx 


/ J sinOdOdx 
jl/2 


1/2 


X 100, 


N,Nb% 


(5.40) 

(5.41) 

(5.42) 


^tns 


xlOO. 


(5.43) 
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with Afij - fi^‘^ - fir fij°^ fij are the exact aad computed solutions respec- 
tively, at the location (i, j) in the numerical grid. and Ng are the total number 
of nodes in the x and 9 directions. As in 1-D, the norms defined above arc also 
used to measure the degree of relaxation of an arbitrary di.stribution at different 
times. In this case represents the Maxwellian to which the initial distribution 
function is expected to evolve. When used for this purpose, the norms are demoted 
by Ln 

5.3.2 Performance of the Quadrilateral Finite Elements 

As in Chapter 3, the performance of the various types of quadrilateral finite elements 
can be evaluated based on the errors introduced by numerical evolution an<l by 
studying the relaxation of a Gaussian distribution function. In the oise of 2-D 
problems, as the order of the element increases, say from linear to quadratic etc, 
the CPU time also increases significantly. A major portion of the time is usually 
spent in solving the banded algebraic system of equations. Nevertheless higher order 
elements produce more accurate results. This fact has already been established 
with 1-D finite elements. 

The size of the banded system also depends upon the way the nodes in the finite 
element mesh are numbered globally, either row wise or column wise. In the present 
problem the semicircular domain is discretized into a finite number of elements. 
Starting from x = 0, the nodes in the 9 direction are numbered first in the asc(;nding 
order from ^ = 0 to tt. This has some advantage as the number of grid points in 
the 9 direction is usually less than that in the x direction. This procedure reduces 
the bandwidth of the algebraic system. The bandwidth is determined from 

iVSW= max [\NODE{IE,i) - NODE{IE,j)\ + ll (5-44) 
1<IE<N^N9 

1< iQ <iVe 

where the matrix NODE{.,.) (known as the connectivity matrix) contains the in- 
formation about the global node numbers in the lEf' element, AT, is the number 
of nodes in each element and ij are the locations of the nodes. Total number of 
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elements {Nx Ng) is the product of the number of elements in the x and 9 directions. 
It can be deduced that the bandwidth is minimum for the linear elements with 4 
nodes and maximum for the cubic elements with 16 nodes. Further the bandwidths 
corresponding to the incomplete elements are less than the corresponding complete 
elements. 

Spatial Convergence 

The relaxation of an arbitrary distribution of the form 

/(x, 6) — a exp(— &(x - a)^ — c(sin 9 — df) (5.45) 

has been studied to estimate the accuracy of the two dimensional quadrilateral finite 
elements. The constant c is varied from 2 to 20 with d = 0. In all the cases it was 
found that the initial distributions eventually reach to the steady state Maxwellian. 
Numerical results for one such case (c = 20) is shown in Fig. 5.2. The relaxation 
measures Tqoi ^2, and ^rm» are shown at 6t (t - relaxation time) for the Crank- 
Nicolson scheme. The number of nodes used in the 9 direction is 31 and it is kept 
constant for all the cases. 

The figure indicates that the linear elements with less number of nodes behave 
differently compared to the quadratic and cubic elements. However as the number 
of nodes is increased, the measures show a trend to move towards those obtained 
using the higher order elements. This behavior coincides with a similar observation 
made earlier in Section 3.4 with 1-D elements. The measures at 6r are different 
corresponding to the different types of elements. Such differences will become much 
less if the number of elements in both x and 9 directions are increased further. It 
is clear that the relaxation norms obtained by the use of incomplete quadratic and 
cubic elements are not much different from those of the corresponding complete 
elements. This is a useful observation as the bandwidth associated with the incom- 
plete elements is less, as observed above. Similar results were observed with the 
fully implicit scheme too. It should be observed from the figures, that while using 
higher order elements, there is little advantage in going from 31 to 91 nodes. The 
degree of relaxation shown is more or less the same in all the four norms. 
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Grid spacing 


1 



Grid spacing 




Figure 5.2: Spatial convergence of various types of quadrilateral finite elements 
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Time step 



Tine step 


Figure 5.3: Temporal convergence of the quadrilateral finite elements 


Temporal Convergence 

Temporal convergence of the five types of quadrilateral finite elements is also studied 
by observing the relaxation of an eurbitrary initial distribution. The results are 
illustrated in Fig. 5.3. The relaxation measures have been obtained with d = 0 in 
eq. (5.45) for a mesh of 61x31 (x, $) nodes. Time step Ar is varied from 1/2 to 
1/32 r. The measures Loo and L 2 , shown in the same figure for the Crank-Nicolson 
scheme, indicate that when the time step is large (Ar = 1/2), all the elements 
slow down the relaxation. As the step size is made small, the linear elements show 
slightly different behavior between At = 1/4 and 1/16. Eventually, for small At, 
they yield the same degree of relaxation as obtained using higher order elements. 

Interestingly the incomplete elements have shown almost the same result as that 
obtained by the use of complete elements. When incomplete elements are used, the 
values at the interior nodes are obtained by interpolation. However, these values 
are needed only for the calculation of Fokker-Planck coefl&cients. Since the latter 
involve only the integrals, the error is relatively small. 

Further the computational time is also an important aspect in judging the per- 
formance of the various elements. The relative CPU times required for 100 time 
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Table 5.1: CPU times for the quadrilateral finite elements with respect to the linear 
elements 


No. 

Element type 

Time 

1 

Linear 

1.00 

2 

Quadratic 

1.99 

3 

Quadratic (I) 

1.09 

4 

Cubic 

3.92 

5 

Cubic (I) 

1.02 


teps on a mesh size of 61x31 nodes are given in Table 5.1. It is natural that the 
;omplete cubic elements, whose bandwidth (of the algebraic system) is maximum 
equire the maximum time. Next in the row in terms of the computation time 
m the complete quadratic elements. However the incomplete elements take ap- 
proximately the same time as that of the linear elements. This is true of both the 
uadratic and the cubic elements. In other words, roughly for the same amoimt of 
omputational effort as required by the linear elements, it is possible to get more 
ccurate results using the incomplete higher order elements. 


lumerical Oscillations 

he problem of the distribution function / becoming negative is investigated with 
simple example. The number of nodes at which / is negative during the evolution 
a highly peaked anisotropic distribution at each time step is shown in Fig. 5.4. 
he results are shown in fractions which are obtained as the ratio of the number of 
Jgative values to the total number of nodes in the mesh. The fractions show clearly 
at when the mesh is made finer from 61x31 (Fig. 5.4a) to 91x31 (Fig. 5.4b), the 
cillations die out faster (in about 3 t in Fig. 5.4b as compared to 5r in Fig. 5.4a). 
^ther fine meshes also reduce the fractional negative values (Chang and Cooper, 
70; McCoy, Mirin and Killeen, 1981; Harrison, 1988). The incomplete elements, 
.reduce somewhat more negative values in the solution which take a little longer 
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a) 61 X 31 mesh b) 91 x 31 mesh 

Figure 5.4: The distribution of fractional negative values of / in a 2-D evolution 


to die out. However it can be visuaUzed that once the spectrum becomes more or 
less smooth after a period of time all the elements show a consistent behavior, with 
incomplete elements yielding almost the same accuracy as the complete elements at 

a lesser cost. 

In 1-D calculations in Chapter 3, fraction of negative / in a multispecies case 
in 250 steps (0.5 s) was presented in Table 3.4. The stepwsie distribution is shown 
in Fig. 5.5 for TV = 61 and 91 nodes. It is observed that coarse grids bring in 
fractionally more negative values in the solution. It is seen that the quadratic 
elements are better than the linear and cubic elements in terms of avoiding the 
spurious oscillations. These oscillations die out after a certain tune. 

5.4 Finite Difference Solutions 

In the finite difference method the computational domain it represented by I and J 
grid points in the * and « directions respectively. The derivatives are approximated 
at each grid point accurate to s^ond order in both directions. The resulting alg^ 
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a) 61 nodes 


b) 91 nodes 


Figure 5.5: The distribution of fractional negative values of / in a 1-D multispecies 
problem 


^raic equations are solved for the distribution functions at these grid points. The 
lize of the system to be solved at every time step is dictated by the choice of the 
ime stepping scheme. 

The 2-D Fokker-Planck equation considered for the finite difference method of 
olution is as given in eq. (2.37). As elsewhere in this work, the dimensionless ve- 
3city and time scales for each species are based on the corresponding characteristic 
alues. The x and 9 derivatives are discretized at a node {i,j) in an unequally 
paced grid using the standard central finite differences (McCoy, Mirin and Killeen, 
981). They are 

» (5.46) 

hj 




9W) 

dx 


(5.47) 
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dx V de, 




1 

2Axi 


Ci+i ,• ( ithtlzJi+hti 

. " I 2A0, 


d{Df) 


de 


hj 


de l^9x. 


- a_i. ( 

V 2A% )\ ’ 


iR. 1 / 

2A% \ m;; 


Fij-1 


de [ de, 


>1? 


Aej 


Fij+l/2 


fi+l,j-l fi-l,j-l 
< 2Axi ^ 


Ae 






Vl/2 


/ /y /»,j-A 

I A0j_1/2 / 


where 


^1/2,2 = 


Ae, 


2 (^2+1 ~ ^ 2 - 1)5 


A^y±l/2 = ±(0^±l-0^), 

Fi,j±i/2 — 2^Fij + Fi^j±i)^ 


and the Axj and Axj+ 1/2 are defined in eqs. (3.48)-(3.49). 
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(5.48) 

(5.49) 


(5.50) 


(5.51) 

(5.52) 

(5.53) 

(5.54) 

(5.55) 


5.4.1 ©-method of Solution in 2-D 

The 2-D Fokker-Planck equations can be integrated in time by different schemes 
including the one step 0-method. As pointed out earlier the time stepping schemes 
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that have already been successfully implemented in the existing codes are, fully 
implicit method, alternating direction implicit method, and implicit operator split- 
ting (McCoy, Mirin and Killeen, 1981; McKenzie, O’Brien and Cox, 1991). In the 
present work, two schemes have been tried, namely the 0 -mcthod and the alter- 
nating direction implicit method. As in 1 -D, the approximation at the boundaries 
can be purely numerical (similar to eq. (3.61)) or based on the analytical limits at 
X = 0, and 0 = 0 and tt derived in Chapter 2 (eqs. (2.121)-(2.125), and (2.133)). 


The 0-Method 

The 0-method has already been introduced in Chapter 3. When this method is 
used, the 2-D equations are solved at all the grid points at the same time. The 
discretization procedure which is second order accurate in space is a nine-point 
difference algorithm. The resulting algebraic equation at each node (i,J) can be 
written as 


T?- Vri-I + = si- (5.66) 

The nine banded structure of the above equation is the consequence of nine-point 
lifference algorithm. The coefficients a ,/?, 7 and 6 are 



-0ar^ 

4-1 = - e ; 

— 

- 0 ci -^ 

ai = 

1 

CD 

/ 3 i =- 05 J + ^; 

li = 

-04, 

II 

+ 

- 0 0^+^ ; 

= - 0 ; 

ti*' = 

- 0 4 ^\ 


(5.57) 


Si = 


(l-0)(arVr+i,;_i+a^/f4 


»+i,j 




Jii 


+IJ+1 
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with 


a: 


al 






j-i 




C^i 


4 


+1 


4 


+1 


+ <4 flu + cj-^Vr-u+i) 


7-“l 


X 


/orn , ^tyJ-1 
^44-1 4 1 


j4AxiA9j \ sinOjJ ’ 


1 fAj^ ^ b.“h,2; 


x?Axi V 2 Axi 


= + 


+ 1 / 2 , 


1 i?*,- 


xj 4 Axi Adj 


, ^+i,i + 


hJ+1 


sin 0j J ’ 


-j. 


A"j-1 ^ •^;-l/2 


xjA9j sin^ y 2 AOj^ij^J ’ 


a:.- 

4 


Axi \Axi+i/2 Axi-i/2j 

1 ( FTjn/2 ^ FTj-y^' 


+ 


+ 


A9j \A9j^y2 A9j^i/2j\ 


1 /'a±i+i5±jzi'i 


X 


j A6j sin0 \ 2 A^j+i/2 J ’ 

1 


X 


% I 


i|4AxiA0j \ sin9j J ' 

1 / ^r-u I 


x^Axi 2 ' Axj_i/2, 

1 /r” +^’ 

^i-lj ■*" 


xj4AxjA% V * sinSj 
The final system is of the form 


Ar+^ = B 


(5.58) 

(5.59) 

(5.60) 

(5.61) 

(5.62) 

(5.63) 

(5.64) 

(5.65) 

(5.66) 

(5.67) 

(5.68) 
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20.0 40.0 60.0 80.0 100.0 20.0 40.0 60.0 80.0 100.0 

Grid points (N) Grid points (N) 


Figure 5.6: Comparisoa of numerical errors using different boundary approxima- 
tions 

where is a vector with I times J components corresponding to grid points 
arranged row wise or column wise and A ia s nine banded matrix. 

'^Jumerical Results 

figure 5.6 presents the errors introduced due to numerical evolution using the two 
lifferent types of approximation at the boundaries, x = 0 and 9 = 0 and tt. The 
esults are presented for different number of grid points. The curves show that 
he approximation based on the analytical limits of the operator at the boundaries 
naproves the Coo error. The ems errors are almost the same in the two cases. This 
an be expected as the approximation at the boundaries effects only a few equations 
1 the system (5.68), Convergence with increasing number of grid points is clearly 
rought out in these figures. 

The relaxation of a Gaussian input distribution function in the form of eq. (5.45) 
as also been studied using both types of approximations at the boundary. The 
rder of the Legendre polynomial used is Z = 4. The relaxation measures Too and L 2 
re given in Table 5.2 for a fixed Ng = 31. The measures show that both numerical 
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Table 5.2: Relaxation norms for a Gaussian input 


N 

Numerical 

•^00 

Analytical 
■^00 I^rms 

31 

1.514 

1.361 

3.288 

3.131 

41 

1.714 

1.436 

2.413 

2.273 

51 

1.854 

1.623 

2.162 

2.009 

61 

1.961 

1.755 

2.081 

1.918 

71 

2.045 

1.845 

2.059 

1.886 

81 

2.110 

1.906 

2.062 

1.876 

91 

2.162 

1.949 

2.074 

1.876 


and analytical approximations yield more or less the same degree of relaxation with 
Ns >71. However when Ns is below 71, the evolution based on the purely numerical 
treatment at the boundaries is faster. Thus for coarser grid there is some ment 
in applying the analytical method at the boundaries. In all these cases the imtial 
distribution approaches the steady state Maxwellian. 


5.4.2 Alternating Direction Implicit Method 

In this method, the equation is advanced in time in two steps. The equation b 
treated as expUcit in d and impUcit in . direction in the first hrff time step 
n to n + 1/2. In the second half time step, that is from n + to n +1, th 
equation is impUcit in . direction and explicit in .. “ethod rs s able ^d 

second order accurate in time. The advantage of this method .s that tt ^ 

tridiagonal systems to advance one step in time. Therefore it is mu* fcter 

e-method. The tridiagon^syst^of equations 

large matrices need not be stored or inverted, as m the case 
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equation for the first and second half time steps can be written as 

9 /^nyn+l/2 ^ 4.(7*^ 




jfj 


.+1 /n-M/2 


•fij JiJ 
At/2 


+ 


cc? sin 9j 




hj 




+ 


a;? sin^j 


d9[ ^ dx de 


(5.69) 


n f In 


(5.70) 

».i 


The mixed derivative terms are treated explicitly. 


'lumerical Results 


t is known that the alternating direction implicit method is stable and has been in 
se to solve the Fokker-Planck equations. The errors caused by the numerical and 
nalytical approximations at the boundaries are shown in Table 5.3 for two species, 
'he spacing in the x direction is varied from 1/30 to 1/90. The errors show that 
nalytical approximation has reduced the Coo errors in both the species. The rms 
rrors are seen to be comparable. This is similar to the results obtained by the 
•-method. 

The results of relaxation of a Gaussian using the alternating direction implicit 
lethod are presented in Table 5.4 for two Ne values (31 and 61). The grid spacing 
. the X direction is varied from 1/30 to 1/90. While there is not a significant 
fference between Ng = 31 and 61, a uniform convergent trend is seen with both 
pes of boundary approximations. 
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Table 5.3: Comparison of errors (percent) with different boundary approximations 
in the alternating direction method for two species (At = 1/16) 


■ 

First species 

Second species 

N 

Numerical 

Analytical 

Numerical 

Analytical 


^00 

^ms 

^00 

^rms 


^ms 

^C» 

^ms 

31 

6.403 

2.263 

3.207 

3.147 

3.771 

1.879 

2.162 

1.957 

41 

2.838 

1.174 

1.402 

1.311 

2.368 

1.054 

1.013 

0.858 

51 

1.611 

0.709 

0.762 

0.686 

1.615 

0.665 

0.560 

0.503 

61 

1.151 

0.468 

0.490 

0.436 

1.157 

0.453 

0.404 

0.356 

71 

0.850 

0.329 

0.358 

0.320 

0.856 

0.328 

0.312 

0.282 

81 

0.641 

0.246 

0.287 

0.261 

0.647 

0.253 

0.261 

0.240 

91 

0.490 

0.198 

0.247 

0.227 

0.495 

0.209 

0.232 

0.215 


Table 5.4: Relaxation norms for a Gaussian distribution with different boundary 
approximations (At = 1/16) 



Theta grid i 

joints = 

31 

Theta grid i 

joints = 

61 

1 

N 

Numerical 

Analytical 

Numerical 

Analytical 1 


LqO 

L'rms 

■^00 

Lrms 

Loo 

Lrms 

Loo 

Lrms 

31 

3.064 

1.500 

2.407 

2.249 

2.973 

1.451 

2.471 

2.325 

41 

1.703 

0.915 

1.470 

1.320 

1.544 

0.839 

1.535 

1.397 

51 

0.973 

0.782 

1.144 

1.027 

0.802 

0.715 

1.211 

1.104 

61 

0.915 

0.791 

1.031 

0.920 

0.856 

0.739 

1.097 

0.998 

71 

0.960 

0.830 

0.981 

0.880 

0.904 

0.787 

1.050 

0.957 

81 

0.997 

0.866 

0.963 

0.865 

0.941 

0.826 

1.033 

0.943 

91 

0.956 

0.828 

0.960 

0.862 

0.891 

0.781 

1.029 

0.937 
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5.5 Background Species and Energy Transfer 

In plasmas fast ions which arise out of neutral beam injection and as charged fusion 
reaction products transfer their energy to the background ions and electrons. In 
fusion related plasmas the self equilibration time of electrons is of the order of 
microseconds. Therefore the electron population thermali'/os rapidly on a faster 
time scale compared to the ions. In general the background species, which are at 
lower temperatures than the high energy species, can be treated as Maxwellian. 
The interaction between fast ions and the cold ions and electrons (Maxwellian) can 
be taken into account easily through the Fokker-Planck coefficients as described 
in Section 2.3.2. This improves the efficiency of computations. The two possible 
approaches are discussed below. 

5.5.1 Fixed Background Species 

Two representative runs are shown in Fig. 5.7 for a single fixed background species. 
?'igure 5.7a is shown for the case of fast ion density, ni 5 (fast) = 0.2 np (background), 
vith an initial fast ion energy of 75 keV. The initial energy of the backgroimd ions 
s 15 keV. The results have been obtained using the 1-D finite element code in the 
ibsence of source and loss terms. Since the energy transferred to the cold species is 
gnored, the fast ions eventually reach the steady state value of 15 keV. A similar 
un with two fast ion species (deuterons and tritons) with deasities n£)(fest) = 
ir(fast) = 0.2 ud (background) is shown in Fig. 5.7b. The transfer of energy from 
he hot species is shown for about 0.5 s. 

.5.2 Energy Transfer between the Species 

he temperature of the background species can be modified if the energy transfer 
ites are known correctly. When arbitrary number of species are modeled the energy 
ansfer rates should be known individually. The rate of energy density gained 
• lost by the species o due to the collision with species b is equal to the second 
oment of the Coulomb operator which represents the collision between species a 
id b (McCoy, Mirin and Killeen, 1981; Killeen, Kerbel, et al, 1986) and is given 
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(a) 


Figure 5.7: Energy loss from a general 
species 



0.0 0.1 0.2 0.3 0.4 0.6 

Time I secs) 

(b) 

in the presence of a fixed badcground 


by 

Qcb = TrmaJ J vUm9dv(W, (5-71) 

= - 2Trmara j sinOdB j G\vdv, (5 • "^2) 

where is the corresponding portion of Ga which depends on the species b. Ga is 
defined in eq. (2.9). The transfer rate in terms of the dimensionless variables is 

sin9d^ G^x^dXa, (5-73) 

where C = nl/va- Va and n„ are the velocity and density normaUzation constants of 
the species a. When both species a and fe are isotropic, the above equation reduces 
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to 


Qab 



i: 


OVa\ 


+ uHio (/‘) 


u 



(5.74) 


Further, if the species b is MfiKwellian, and the species a is isotropic, then 


Qofr = 





% dfa\ 
rribvl dx j 



(5.75) 


where the functionals M,N, and R are as defined in eqs. (2.56)-(2.59). Tj, is the 
temperature of the species b. It may be noted that when many species are considered 
in the simulation, calculation of energy transfer rates at each time step would require 
:onsiderable run time. Provisions have been made in the computer programs to 
:alculate the transfer rates. If only one background species is considered, its energy 
:an be updated easily without the knowledge of the energy transfer rates. Since the 
inergy conservation can be maintained, both in the presence and absence of source 
ind loss terms, through the procedure outlined in Section 3.5.1, the difference in the 
otal energy of the general species can be added to the background species. This 
las been done particularly to treat electrons as background species, if desired. 


Fig. 5.8 illustrates the energy transfer from a general fast ion species to a back- 
round species where the energy is updated as explained above. Fast deuterons 
dth a density of nn = 0.2 ny sharing their energy with backgroimd tritons is shown 
1 Fig. 5.8a. It can be seen that the energies of both the species approach the 
verage value of 25 keV in about 0.5 s. Another typical situation is illustrated in 
ig. 5.8b. In this case two general species, deuterons and tritons, are allowed to 
ain energy from the background deuterons. These results are shown to establish 
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(a) 


(b) 


Figure 5.8: Energy transfer between general species and a varying background 
species 


that the models used for the background species adequately and accurately account 
for the energy exchange. 


To study the possibilities of treating electrons as background species, a deuteron- 
electron plasma was considered with an average energy of 45 keV. Initially the fast 
ions having 75 keV were allowed to relax their energy with background electrons 
(15 keV). The energy exchange is shown in Fig. 5.9. The deuterons and electrons 
are shown in solid and broken lines. The same study was performed treating both 
electrons and ions as general species. The results are shown as circles and triangles 
in the same figure. It can be seen that treating electrons as background species 
may not raise any serious errors in the numerical simulations. This approach will 
reduce the computation time considerably because the Fokker-Planck equations can 
be integrated on the ion time scale throughout the evolution. 





124 


CHAPTERS, numerical SOLUTION OF 2-D EQUATIONS 



0.0 0.5 1.0 1.5 2.0 2.5 

Time (secs) 

Figure 5.9; Energy transfer between ions and electrons 

5.6 Summary 

Fokker-Planck equations in 2-D have been solved by the finite element and the 
inite difference methods. The finite element computational domain is a semi circle 
md is discretized into a number quadrilateral elements. The equations have been 
ntegrated in time by the one step 0 method. The relative performance of the 
various finite elements has been assessed. Alternating direction implicit method has 
ilso been applied in the finite difference method. The main points of the present 
hapter can be summetrized as follows: 

i. Multi-velocity scale Fokker-Planck equation in two dimensions have been 
solved using five types of quadrilateral elements namely, linear, and complete 
and incomplete quadratic and cubic finite elements. 

ii. Based on the L„ relaxation measures for Gaussian distribution ftmetions, it 
is found that the performance of the quadratic and cubic finite elements is 
better than the linear elements. Further, the incomplete elements are found 
to compare well with the corresponding complete elements. 
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iii. The quadratic elements are preferable among the higher order elements. The 
problem of numerical oscillations in the case of quadratic elements is much less 
compared to the cubic elements. For approximately the same computational 
effort as that for the linear elements, incomplete elements can be used to solve 
the equations to obtain more accurate results. 

iv. The analytical limi ts of the collision operator at the botmdaries, x = 0, and 
9 = 0 and tt, have been implemented for the 2-D Fokker-Planck equations. 

V. Background species with fixed as well as varying temperature have been mod- 
eled in the finite element method. Numerical experiments show that the 
electrons can generally be treated as background species unless the electron 
distribution is explicitly required. 



Chapter 6 


Fusion Reactivities and Fast Ion 

Interactions 


6.1 Introduction 

As stressed earlier accurate description of plasma behavior requires considerable 
effort in terms of understanding and modeling. A realistic representation of plasma 
within the frame work of kinetic models should take into account the various mecha- 
nisms that occur in plasmas. Therefore, in addition to the modeling of Coulomb in- 
teractions in velocity space using the Fokker-Planck collision operator, other mecha- 
nisms /processes should also be included in the simulation. Even though an accurate 
modeling of the relevant processes may be quite complex, in general representative 
models can always be used in a broad manner in order to ascertain their impact oh 
the distribution functions of charged species in fusion plasmas. 

Shape of the distribution functions of a species is important under non-Maxwellian 
conditions. Changes in the fusion reactivities due to tail formation, and energy 
gain or loss are some of the examples which are affected by the actual shape of the 
distribution functions of the interacting species. Tails are mostly formed due to 
auxiliary heating and the reaction products. Under such scenarios, the numerical 
models that describe these processes should be explored fuUy. This may be achieved 
by studying/varying the different parameters pertaining to each of the source/sink 
models considered in the simulation of fusion plasmas. 
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L28 chapter 6. REACTIVITIES AND FAST IONS 

In this chapter, an improved fusion reaction loss model is presented. Fast ion 
iources, and particle and energy sinks are also included. Fusion reac tivities are 
:alculated in the presence of sources and sinks. The shapes of the distribution 
unctions in some of the transient scenarios are also illustrated. 


>.2 Reaction Losses 


L general form of the Fokker-Planck equation with source and sinks for each species 
an be written as 


dfa 

at 



+ SaF + SaR — 
c 





( 6 . 1 ) 


rhere the first term in the right hand side of the above equation is the Fokker- 
lanck collision term, the second and third terms correspond to fast ion and reaction 
jurces, and the last three terms represent the particle, energy, and reaction sinks, 
or each species a, the eq. (6.1) is solved with the relevant source and loss terms. 


.2.1 Killeen’s Model 

he fuel ions in a multispecies plasma deplete due to fusion reactions. The reaction 
nks are calculated using the loss coefficient. The loss coefficient for species a due 
I reaction with species b is denoted by Lab- According to Killeen, Kerbel, et al. 
986), Lab is given by 

Tab ~ fbb ^ CV ^abj ( 6 . 2 ) 

iere Ub is the density of the reacting species b, and < erv >a6 is the fusion reactivity 
itween species a and b. The above approach is referred as Killeen’s model in this 
)rk. The reactivity < av > Is defined as 

I M'i'Mu) u. 


(6.3) 
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whore fa &nd /j, are the distribution functions, (t{u) is the reaction cross section 
which is a function of the relative velocity u. The above integral can be evaluated 
easily using the procedures suggested by Marx, Mirin, et al. (1976), and Gordey, 
Marx, ct al. (1978). 

The reaction loss coefficient in eq. (6.2) is independent of energy. This implies 
that the fraction of particles lost of species a is taken to be constant for all Xa (the 
dimensionless velocity scale for species a). 


6,2.2 Improved Model 

In fusion plasmas high energy fuel ions react faster with each other due to higher 
fusion cross-section. The ions in the high energy region are bimit more due to fusion 
reactions than the low energy ions. Therefore the losses at all the energy levels 
need not be uniform. In order to represent the reaction losses more accurately the 
model given in eq. (6.2) can be improved. An improved model is presented here for 
this purpose and also implemented in the codes developed in the present work. The 
reaction losses for the 2-D problems are described in the nondimensional variables 
by using the loss coefficient 

roo fir f^ir f 

Lab{xa,ff) I dxhxl t d9i2 sm9i2 d<j>i 2 fb{xb,^)cr{u)u, (6.4) 

Jo Jo JO 


and in case of 1-D models 

Lab{Xa) = 47rn6 r dxb xl f ddnfbixb) <riu) u e.a (6.5) 

Jo JO 

where the symbols have their usual meaning. The relative velocity u is written as 

U = [vlxl + vfxl - 2 Va Vb Xa Xb COS ^12] (^•®) 

and 012 is the polar angle between the velocity vectors representing Xa and Xb. Here 
Lab is a function of dimensionless velocity Xa and the integrals on the right hand 
sides of eqs. (6.4)-(6.5) are evaluated for each value of a;„. For 2-D problems, it also 

becomes a function of the pitch angle 9. 
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(a) D-T plasma (b) D-He^ plasma 

Figure 6.1: Reaction loss coefficients of fuel ions in D-T and D-He^ plasmas 

The loss coefficients Lab are shown in Fig. 6.1 for D-T (15 keV) and D-JTe^ 
>0 keV) plasmeis. Densities of the species are 10^ m~^ {ud = nr = n). The 
responding n < av > values, as obtained by the Killeen’s model (6.2), are also 
lown in the figure. It is evident from the figures that there will be a significant error 
the loss coefficient is tadcen to be constant over the velocity (shown as n< av > in 
le figure). The model presented here correctly accounts for the fraction of particle 
st at all energies, which is significantly higher at higher energies, as expected. The 
ss coefficient shown in Fig. 6.1 was obtained for Maxwellian distributions. 

In order to see the effect of the shape of the distribution function on the loss 
efficient, distributions with tails are chosen. A typical distribution function with 
tail is shown in Fig. 6.2. The tail has been introduced artificially. The shape 
s been modified to appear similar to the one observed elsewhere (Cordey, 1975; 
ckerton, 1979). The increase in the average (kinetic) temperature due to the tail 
[ess than 10%. The corresponding loss coefficients Lab are shown in Fig. 6.3 for D- 
and D-jETe^ fuel ions. The presence of the tail, while it enhances the reaction rates 
the high energy end, does not drastically effect the fractional losses. However 
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Figure 6.2: Distribution function with a tail 


some effect on the fraction of particles lost in the low energy region is clearly visible 
in the figures. This is due to the interaction between the tail and the slow ions, 
leading to enhanced burn up slow ions. 

In view of the above observations, the reaction loss term included in the Fokker- 
Planck equation (6.1) is of the form 



where the sum over b is for all the relevant reactions, and is the characteristic 
time constant for ions. 

The impact of the improved reaction loss model on the reactivities is studied with 
a D-T plasma. For different temperatures, the deuterons and tritons were allowed 
to react without any other source/sink, over a period of 5.0 s. The results are 
shown in Table 6.1. The reactivities have been obtained (in units of using 

both the Killeen’s model and the model presented above. The differences in the 
reactivities between the two models are presented in terms of a parameter e defined 
as e = (< (TO >K - < <^« >/)/ < >/• Here the subscripts K and I correspond to 
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Table 6.1: Comparison of fusion reactivities in a D-T plasma as obtained by Killeen’s 
and Improved reaction loss models 



f = 1 s 

t = 5 s 


< erv >K 

< av >i 

e in % 

< av >K 

< av >j 

e in % 

5 

0.1282e-22 

0.1276e-22 

0.05 

0.9567e-23 

0.9370e-23 

2.1 

10 

0.1120e-21 

0.1096e-21 

2.19 

0.1032e-21 

0.9288e-22 

11.1 

15 

0.2735e-21 

0.2658e-21 

2.90 

0.2649e-21 

0.2302e-21 

15.1 

20 

0.4346e-21 

0.4234e-21 

2.60 

0.4281e-21 

0.3771e-21 

13.5 

25 

0.5689e-21 

0.5576e-21 

2.03 

0.5643e-21 

0.5121e-21 

10.2 

30 

0.6723e-21 

0.6630e-21 

1.40 

0.6693e-21 

0.6255e-21 

7.0 


where S denotes the birth rate of the corresponding fast ions. The distributions of 
the source ions are denoted by wj, and tej respectively, and are taken as a highly 
peaked Gaussian functions, with peaks corresponding to the birth energies of the 
fast ions. 

To study the change in the distribution functions of interacting species in a 
plasma, the evolution of a deuteron-electron plasma with fast ions (source) is fol- 
lowed for 2 s. This is done by treating the fast ions as a general species but with 
zero initial density. The ion source is added to this species. The electron popula- 
tion is modified to satisfy the neutrality condition. The electron density is adjusted 
suitably without changing the average energy of electrons at each time step. The 
distribution functions are shown in Fig. 6.4 for the fast ions and electrons.^ The 
results are given for the fast ion energy of 200 keV and birth rate 0.5 x lO^" m ^ s . 
The Maxwellian distributions corresponding to the temperatures at 2 s are also 
shown in the same figure by dashed curves. It can be observed (Fig. 6.4a) that 
sUghtly more number of particles are in the low energy region compared to that 
of Maxwellian, and a slight depletion is at the tail end. If the fast ion source is 
switched off after 2 s, the Coulomb collisions will eventually bring the species into 
equilibrium. It may be observed that in general the shapes of the distnbution 
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(a) 


(b) 


Figure 6.4: Distributions of fast ions and electrons in a D-e plasma 


iction remain close to Maxwellian after sufficient time. 

Figure 6.4b shows the distribution function for electrons. It is evident from the 
ire that the distribution is not any different from the Maxwellian, which is also 
iwn in the same figure. This is because the electrons interact and thermalize on 
ister time scale compared to the ions. Therefore electrons can as well be treated 
a background species. To appreciate this further, the above experiment was 
eated for different energies (Bp)of fast ions by treating the electrons as a general 
I a background species (see Section 5.5). The results are shown in Table 6.2. 
an be seen that both the treatments have shown almost the same results. The 
:ribution function also confirms this (Fig. 6.4b). The net fraction of energy 
isferred from fast ions to electrons is in the range of 30-35%. These observations 
gest that for many applications electrons can be treated as background species, 
with varying energy. 
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Table 6.2: Comparison of energies of different species obtained by treating electrons 
as a general and a background species (Gen-General and Bgd-Background). 


t 

Ef 

Deuterons 

Fast 

ions 

Electrons 

in s 

in keV 

Gen 

Bgd 

Gen 

Bgd 

Gen 

Bgd 


50 

14.50 

14.15 

14.98 

14.73 

12.46 

12.96 

1.0 

100 

24.92 

24.38 

27.63 

27.42 

19.05 

19.79 


150 

35.38 

34.66 

42.69 

42.58 

24.64 

25.59 


200 

45.56 

44.62 

60.22 

60.21 

29.48 

30.64 


50 

19.55 

19.41 

19.65 

19.60 

17.74 

18.06 

2.0 

100 

35.60 

35.81 

36.55 

36.95 

29.37 

29.40 


150 

52.27 

52.96 

55.03 

56.02 

39.47 

39.02 


200 

69.19 

70.46 

74.94 

76.56 

48.39 

47.39 


6.3.2 Reaction Products 

The generation of high energy reaction products is also modeled m the same manner 
as described in the above section. In a D-T plasma, the D-T and D-D reactions 
produce respectively 3.52 MeV alphas, and 1.01 MeV tritons and 0.82 MeV He^ 
ioas as reaction products. The generation rates (birth rates) of alphas and tritons 

arc given by 


For alphas Sa — n^nT < >j)t > 
For tritons St — O.Sn^ <(tv >ddi 


(6.9) 

( 6 . 10 ) 


where n denotes the density of respective species. These rates are need m e 
reaction source term (S.n in eq. (6.1)). The distribution functions and for 
the reaction products axe also modeled as highly peaked Gansslans, as m the case 

of fast ions. 

A simulation experiment was conducted to study the evolution of plasma ^ 
prising four general species which include electrons also. The other speaes are 
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sutcrons, tritons, and alphas. All tho spccios aro roproscntcd on thoir c.haractor- 
tic velocity scales. Only reaction source for alphas and the loss of fuel ions due 
) reactions are included in the simulation. Initial density of alphas is taken as 
The transient behaviors of deuterons and alphas at 1.5 s are shown in Fig. 6.5. 
can be seen that the distribution function of deuterons, denoted by solid line 
Fig 6.5a, remains close to Maxwellian. The Maxwellian distribution at the cor- 
sponding temperature is shown in the same figure by dashed curve. It may be 
)ted that, compared to the Maxwellian, the population in the high energy region 
slightly more indicating a tail like formation. This confirms the fact that high 
lergy reaction products do modify the distribution function in the high energy 
gion, although marginally in this case. The distributions of tritons and electrons 
e also found to be close to Maxwellian similar to those shown in Fig. 6.5a. 


The distribution of alphas at 1.5 s (Fig. 6.5b) shows that beyond a certain 
iocity range (v = O.SiJa), the distribution function falls-off rapidly. Thus, it is 
equate to choose just two times the velocity corresponding to the birth energy 
3.52 MeV. The figure also indicates that although alphas are born relatively at a 
y high energy, due to Coulomb interactions some alphas still gain more energy, 
is is mainly due to a — a interactions. Due to transfer of energy from alphas to 
ler species, the buildup of alphas at low energy is evident from the distribution 
iction. 


The change in the fusion reactivity < crti > at different times due to alpha 
ting in the above simulation is shown in Table 6.3. There is a steady increase 
he fiision reactivity up to 2 s. The biurn up is also shown in the same table for 
initial density of njy = 10^ m (njj = nj). The burn up has gone up to 5% in 
due to alpha heating. It can be seen that the fusion reactivity shows decreasing 
id at 2 s. This is due to the fact that fuel ions have gained sufficiently high 
■Sy by this time and further gain in energy does not lead to increase in reactivity 
to the nature of D-T cross-sections. In this simulation, the particle and energy 
mechanisms were not included. 
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(a) 


(b) 


Figure 6.5: Distribution functions at 1.5 s in the presence of reaction products 


Table 6.3: Increase in the fusion reactivity due to alpha heating in a D-T plasma 


time 

s 

< (TV > 

rr? s~^ 

Burn up 

% 

0.171 

0.140e-21 

0.21 

0.266 

0.160e-21 

0.31 

0.485 

0.220e-21 

0.76 

0.775 

0.297e-21 

1.54 

0.972 

0.332e-21 

2.08 

1.529 

0.358e-21 

3.79 

1.937 

0.413e-21 

5.19 

2.000 

0.385e-21 

5.41 
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4 Peirticle Loss 


e particle loss term used in eq. (6.1) is of the form (McCoy, Mirin and Killeen, 
11; Matsuura, Nakao and Kudo, 1992), 




( 6 . 11 ) 


ere Tp(x) is the particle confinement time and is velocity dependent. It is given 
;erms of the dimensionless velocity variable x as 


Tp(a;) = r^max[l,cx]^. (6.12) 

i velocity dependence of the confinement time is achieved by varying the param- 
• p. Tq can be adjusted to keep the total loss of particles constant (Matsuura, 
lao and Kudo, 1992; Matsuura, Nakao, et al, 1993). The constant c is related 
he ratio of the maximum cut-off to the thermal velocity. 

For a single species, when confinement time is kept constant, i.e, j3 = 0, and 
1 only the fast ion source and particle loss terms, the Fokker-Planck equation 
) can be integrated to obtain both density and energy. That is, the density at 
time t is given by 

n(t) = [n(0) - rp Sp] exp-*^^ +Tp Sp, (6. 13) 

re n(0) and Sp are the initial density and the birth rate of the faust ions respec- 
y. The average energy can be calculated by using 

Eit) = ( [n£’(0) - Tp Sp Ep] exp-*/^^ +rp SpEp)/n{t), (6. 14) 

re the values at i = 0 are specified. 

. he effect of parameter is shown in Fig. 6.6. Computational results for a 
e species are shown for /3 = 0 to 4. The fast ions are added to the bulk ions at 
eoi S = 0.5 X 10^ m ^ s The results are shown for Tq = 2.0 s. The straight 
corresponding to /? = 0 in Fig 6.6a shows that the particle loss is balanced by 
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(a) (b) 

Figure 6.6: Variation in density and energy for a single species with particle sink 


the fast ion source, as is evident from eq. (6.13). As j3 increases, the number of 
particles lost in the high energy region becomes less compared to the leakage of the 
low energy particles. This results in a net increase in density as shown in the figttre. 
The average energy of the species is shown in Fig. 6.6b, and shows a relatively less 
variation with the parameter /3. 

Single species calculations were also made to observe the shape of the distribu- 
tion functions during evolution in the presence of particle sinks. A fast ion source 
with energy Ff = 200 keV is added to the bulk deuterons. The results are shown 
in Fig 6.7 for two cases: (i) fast ion source only (Fig. 6.7a), and (ii) fast ion source 
and particle sink /3 = 1 (Fig. 6.7b). The corresponding Maxwellians are also shown 
in the same figures by dashed curves. Reactivity calculations performed for a D-T 
plasma at different temperatures, under the above mentioned conditions, have not 
shown appreciable differences from those of the corresponding MaxwelUan values. 
Thus if the energy loss term is not included, the effect of the parameter /3 on the 
distribution function or the reactivities appears to be small. 
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Figxrre 6.7: Distribution function of deuterons in the presence of particle sink 


5.5 Energy Loss 


he energy loss term used in the dimensionless velocity variable is given by (Mat- 
uura, Nakao and Kudo, 1992) 



ii.[ f] 

x^dx [ 2 t £( x ) ’ 


(6.15) 


here T£ represents the velocity dependent 


energy confinement time and is written 


Te(®) -- Tq max[l, cx]^. (6.16) 

can be adjusted to keep the velocity integrated energy loss constant. 

ugle species, the Fokker-Planck equation with fast ion source and particle 
energy sinks can be integrated to obtain energy and density for /3 = 0. 

When fast ion source and only energy sink terms are present: 
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E{t) = {[nE{0) - teSfEf] + tfSfEf) /n(t). (6.17) 


Here n(t) is the density at time t as obtained in eq. (6.13). 
ii. Wiien fast ion source, and both particle and energy sinks are included: 


E{t) = {[nEiO) - t.jjSfEf] 4- r,ffSFEF)/nit), (6.18) 


where n(t) is as in eq. (6.13), and Te// is given by 


_!_ = jL + JL. (6.19) 

Te// Tp T£ 

The effect of the parameter ^ on the energy loss in a single species in the presence 
of fast ion source is shown in Fig. 6.8 for rp = T£ = 2.0 s. It can be seen that for 
/? > 2, there is not much change in the average energy of the species. This also 
does not change the distribution functions appreciably. This fact can be seen from 
Table 6.4 in which the relative increase in the fusion reactivities for different /3’s 
and also for different EpS, the fast ion energies, are given. The relative < av > 
values are shown in terms of < an > for ^ = 0. These values are shown at 1 s 
and 2 s for two different fast ion energies. It is observed that the enhancement in 
the reactivity is about 18% for /3 = 0 to 2, compared to 1% for ^ = 3 to 4. The 
comparisons of these reactivities with the respective MaxwelUan values at the same 
kinetic temperature indicate negligible differences. This indicates that the fast ions, 
added to the bulk ions, obtain a distribution close to MaxwelUan in 1 to 2 s. 


6.6 Summary 

An improved reaction loss model which correctly calculates the particle loss due 
to fusion reactions at different energies has been implemented. Its impact on reac- 
tivity at different temperature is shown. The sources and sinks in the srmulatron 
studies are fast ion sources including the reaction products, and particle and ener^ 
losses. It is found that the high energy fusion product alphas can be represen e 
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Figure 6.8: Vaxiation in energy of the species (bulk+fast) in the presence of sinks 


able 6.4: Effect of parameter /? on < cru > in time dependent calculations of D-T 
lasmas with particle and energy sinWg 


p 

Ep = 

100 keV 

Ef = 

150 keV 


1 s 

2 s 

1 s 

2s 

0.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.144 

1.148 

1.078 

1.050 

2.0 

1.185 

1.181 

1.091 

1.054 

3.0 

1.194 

1.183 

1.092 

1.054 

4.0 

1.193 

1.180 

1.091 

1.053 
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a(l<KiUfit('ly in the velocity space by choosing a cut-off twice the velocity at their 
birth eiK'rgy. The distribution function for alphas indicates that a — a interactions 
yield some a particles above their birth energy. The simulation studies conducted 
by treating electrons as a general species also suggest that electrons can be treated 
as a background species with varying energy. The shapes of the distribution func- 
tions and their role in the calculation of fusion reactivities in the presence of velocity 
d<'i)endent loss mechanisms are studied. Generally, the difference in the bulk ion 
r(^activities compared to the corresponding Maxwellian reactivities at the same tem- 
p<'ratur<'s is not significant. 



Chapter 7 


Concluding Remarks 


7.1 Summary and Conclusions 

In tho prcst'.nt work various computational aspects pertaining to the numerical so- 
lutioas of Fokker-Planck equations have been investigated. The central idea behind 
this effort is to improve the method of modeling ^ multispecies fusion plasmas, 
and also to enhance the accuracy in the computed solutions, where possible. The 
contributioas of the present work are detailed in the following paragraphs: 


The concept of multi- velocity and time scales in Fokker-Planck equations us- 
ing characteristic velocities and relaxation times for each species has been 
presented. The dimensionless Fokker-Planck equations and the Rosenbluth 
potentials in 1-D and 2-D velocity space are given. These equations can be 
used for an arbitrary number of general and background species. In the case of 
background species, separate expressions for the coefficients have been arrived 
at. These equations are solved by the finite element and the finite difference 
methods, using Crank-Nicolson and fully implicit time stepping schemes. 

One dimensional multi-velocity scale Fokker-Planck equations have been solved 
by the Galerkiu finite element method using Unear, quadratic, and mbic e e- 
ments. It is found that the quadratic aud cubic elements perform better m 
Unear elements in terms of computation^ efficiency and convergence vnth 
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mesh refinement and time increment. In the case of 2-D, the equations are 
solved using five types quadrilateral elements namely, linear, and complete 
and incomplete quadratic and cubic elements. Based on the relaxation norms 
for Gaussian distribution functions, it is again concluded that the performance 
of the quadratic and cubic elements is better as compared to the linear ele- 
ments. Further, the incomplete elements are found to compare well with the 
corresponding complete elements. It is seen that for approximately the same 
computational effort as that for the linear elements, incomplete elements can 
be used to solve the equations to obtain more accurate results. 

• The distribution fxmctions becoming negative (erroneous numerical oscilla- 
tions) at the tail end has been studied in the multi-velocity scale finite element 
solutions. It is observed that the choices of the characteristic velocity nor- 
malization constants (cut-offs) reduce the numerical oscillations significantly. 
This fact is more evident from the numerical experiments in which consider- 
able reduction in the negative values was seen when the velocity cut-offs were 
refined/changed during evolution. 

• Electrons have been treated as a general species together with ions. In the 
multi-velocity scale model, the energy transfer between ions and electrons 
can be studied using the same number of grid points. The results of Spitzer 
model which is valid under Maxwellian approximation are found to be in good 
agreement with the corresponding Fokker-Planck solutions. It may be added 
that the finite element results show better agreement with the Spitzer model 
compared to the finite difference method. 

• Background species with fixed 2 is well as varying temperature have been mod- 
eled in the finite element method. Numerical results indicate that electrons 
can be treated as background species unless the electrons distribution func- 
tion is explicitly required. It is observed in the simulation experiments that 
the transient distribution fimction of electrons remains close to Maxwellian in 
general. 
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• The Fokker-Planck collision operators both in 1-D and 2-D have singularities 
at the boundaries. The singulairities have been removed analytically and the 
non-singuleir operators to be applied at the boundaries have been given. These 
operators have been implemented in both 1-D and 2-D. In case of 2-D it 
has been implemented both in the © -method and the alternating direction 
implicit method. 

• Comparison of the efficiency of the finite element method with the finite dif- 
ference method indicates that for the same computational effort, the former 
method performs better than the latter. The experience of using electrons as 
a general species also supports this as seen in the experiments with Spitzer 
model. 

• In the multi- velocity scale approach the nodes for different species do not cor- 
respond to the same physical velocities and interpolations have to be used in 
the calculation of Rosenbluth potentials for the Fokker-Planck coefficients. In 
order to improve the accuracy, higher order schemes have been used both for 
interpolation and quadrature. It is found that the cubic polynomials based 
Hermite and Lagrangian schemes improve the accuracy in the numerical solu- 
tions as compared to the trapezoidal method. The projection errors are seen 
to be negligible with cubic polynomials. 

• In the present work both densities and total energy are conserved exactly 
during the evolution. This has been achieved by imposing the appropriate 
constraints in density and energy after each time step or a fixed number of 
time steps. This avoids the ciunulative build up of density and energy errors 
in the computed solutions. 

• Local discretization errors have been computed at each node in the compu- 
tational domain. It is shown that by suitably relocating a fixed niimber of 
nodes, the Tnavinrinm error can be reduced by an order of magnitude or even 
more. The reliability of the error estimates used in the grid adaptation is 
checked against the exact errors where possible. 
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• A scheme has been presented to accelerate the time marching in the simu- 
lation. This scheme is based on the energy error in the computed solution. 
The results indicate that the number of time steps, and consequently the 
computational time, can be reduced significantly. 

• An improved reaction loss model which correctly calculates the fractional 
particle loss due to fusion reactions at different energies has been presented. 
It is seen that there can be appreciable error if the loss coefficient is taken to 
be a constant as it is done elsewhere. In addition, fast ions sources, reaction 
products, and sinks represented by velocity dependent particle and energy 
losses are also implemented in the codes developed. 

7.2 Observations and Suggestions 

A few observations and suggestions which may be useful in extending the present 
work are given below: 

• In the present work, in solving the 2-D Fokker-Planck equations using the fi- 
nite element method, the 0-method has been used. The alternating direction 
implicit scheme may have some advantages from the viewpoint of computa- 
tional efficiency and should be tried in the finite element method. 

• Moving finite elements in 2-D should be attempted. An infinite element to 
represent the velocity space beyond the cut-off(s) may also be useful and 
should be tried. 

• In the present work the importance of the convective and diffusive terms in the 
Fokker-Planck equations has not been studied. This may have some impact 
on the numerical solutions and it should be studied in detail. 

• The Fokker-Planck equations using the velocity variables have been employed 
in this work. While some work has been reported in the literature using energy 
as the independent variable, more work is needed to investigate the relative 
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merits of using the velocity or the energy variable in a multispecies plasma 
from a computation viewpoint. 

• RP heating should be included to study the combined auxiliary heating effec- 
tively. Additional source/loss mechanisms/models can be added where needed 
to improve the acciuacy and to make the simulation more realistic. 

• The multi-velocity Fokker-Planck equations should be coupled with radial 
transport models to be able to study the effect of radial transport on heating 
and current drive in magnetically confined plasmas. 
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Appendix A 

Legendre Polynomials and their 

Derivatives 


The Legendre polynomials of order up to 8 and their derivatives are listed here. 
Legendre polynomial of order Z, and its nth derivative are denoted by Pi{6) and 
Pi {9) respectively. 

• 1 = 0 


PoW = 1-0 

P^{9) = cos e 

pi (0) = - sin 9\ Pl{9) = - cos 9\ Pf (6) = sin 9 


• 1 = 2 


P^{9) = (3cos20 + l)/4 
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Pl{9) = (-6sin20)/4 
p2(0) = -3 cos 20 

P|(0) = 6 sin 20 

• 1 = 3 

P^{ 6 ) = (5 cos 30 + 3 cos 0)/8 
PI (0) = (— 15 sin 30 — 3 sin 0)/8 
Pl{e) = - J(45cos30 + 3cos0)/8 

O 

J^(0) = (135 sin 30 + 3 sin 0)/8 

• 1 = 4 

Pji9) = (35 cos 40 + 20 cos 20 + 8) /64 
PM = (-14Osin40-4O sin20)/64 
Pi (^) = (-560 cos 40 - 80 cos 20) /64 
Pl{e) = (2240 sin40 + 16O sin20)/64 

• 1 = 5 

P^i^) = (63 cos 50 + 35 cos 30 + 30 cos 0)/128 
P 5 (^) = (~ 135 sin 50 - 105 sin 30 - 30 sin 0) /128 
Pli^) = (-1575 cos 50 - 315 cos 30 - 30 cos 0) /128 

i^( 0 ) = (7875 sin 50 + 945 sin 30 + 30 sin 0)/128 
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• 1 = 6 

P^{9) = (231 cos 69 + 126 cos 40 + 105 cos 20 + 50) /512 
Pgi(0) = (-1386 sm60- 504 sin40- 210 sin20)/512 
p^{e) = (-8316 cos 60 - 2016 cos 40 - 420 cos 20)/512 
p3(0) = (49896 sin 60 + 8064 sin 40 + 840 sin 20) /512 

• 1 = 7 

^7(0) = (429 cos 70 + 231 cos 50 + 189 cos 30 + 175 cos 0)/lO24 
p^{9) = (-3003 sin 70 - 1155 sin50 - 567 sin30 - 175sin0)/lO24 
p2(0) = (-21021 cos 70 - 5775 cos 50 - 1701 cos 30 - 175 cos 0)/lO24 
p3(0) = (147147 sin 70 + 28875 sin50 + 5103 sin 30 + 175 sin0)/lO24 

• 1 = 8 

Pg (0) = (6435 cos 80 + 3432 cos 60 + 2772 cos 40 + 2520 cos 20 + 1225) 
p^(e) = (-51480 sin 80 - 20592 sin 60 - 11088 sin 40 - 5040 sin 20)/ c 
p2 (0) = (-411840 cos 80 - 123552 cos 60 - 44352 cos 40 
— 10080 cos 20) /c 

if (0) = (3294720 sin 80 + 741312 sin 60 + 177408 sin 40 
+ 20160 sin 20) /c 


where c = 16384. 
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Finite Element Shape Functions 

l-D Shape Functions 

Shape functions for various l-D elements are given below (Dhatt and Touzot, 1985): 

• Linear Elements (2 nodes) 

<i>i = ( 1 - 0/2 
< j >2 = ( 1 - 1 - 0/2 

• Quadratic Elements (3 nodes) 

<f>^ = - 01 - 0/2 
h = (1-0) 

•f>z = 01+0/2 

• Cubic Elements (4 nodes) 

= -(l-0(l-90)/16 

h = 9(l-0)(l-3e)/16 
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cj>, = 9(l-^')(H-3e)/16 
= -(1+0(1 -9^/16 

2-D Shape Functions 

The shape functions ^’s for all the five t)^es of quadrilateral elements considered 
in the present work are listed here (Dhatt and Touzot, 1985; Reddy, 1984). The 
nodes in each element are numbered as shown in Fig. (5.1). The shape function <j)i 
corresponds to the node i in each type of the finite element. 

• Linear Elements (4 nodes) 

= ( 1 - 0 ( 1 - 0/4 
<f>2 = ( 1 + 0 ( 1 - 0/4 
h = (1 + 0(1 + 0/4 
04 = (1-0(1 + 0/4 


• Complete Quadratic Elements (9 nodes) 

01 = (i-0(i-0^V4 

02 = -( 1 - 0 )( 1 - 0 ^ 7/2 
<t>z = -(l+0(l-0^^/4 

0 , = ( l +^)( l _ ^ 2)^/2 

05 = (1+0(1 + 0^W4 

= (l-0)(l+0^/2 



h = -( 1 - 0 ( 1 + »?)^»?/4 

<^9 = (i-OKi-r?^) 

Incomplete Quadratic Elements (8 nodes) 

<l>\ = (1 - 0 ( 1 - 0 ( 1 + ^ + 0/4 
h = (l-0)(l-0/2 
h = -(i + 0(i-0(i-^ + 0/4 

<l>4: = (l+0(l-0)/2 

<^5 = -(1 +0(1+ 0(1 -e- 0/4 

<l>s = (l-0)(l + 0/2 

(j)^ = -(1 -0(1 + 0(1 +^- 0/4 

</>s - (l-0(l-0)/2 
Complete Cubic Elements (16 nodes) 

h = N^iON^iv) 
h = N^iONiiv) 

<j>z = N^iONM 
<j>, = N,{i)N,{ri) 

<^5 = iV4(0^2(0 

<^6 = ^4(0 ^3(0 
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h = N.iON^iv) 

= NsiON^iv) 

h = N2{i)N^{ri) 

^10 = N,{i)N^{ri) 

= N,{i)N^{ri) 

<i>n = N,{^)N2{ri) 

<^13 = N^{^)N2{ri) 

^14 = N,(0N2{v) 

9^15 = NziON^in) 

4>ie = N2(^)N2(v) 

where 

iVi(e) = -(i_^)(1_9^2)/;^6 

N^iO = 9(l-e')(l-3^)/16 
iVaCO = 9(1 -e^)(l + 30/16 
N^iO = -(l+0(l-90)/16 

• Incomplete Cubic Elements (12 nodes) 


9^1 = 9(1-0(1-»7)A/32 

h = 9(1 - 30(1 -0)(1- 0/32 
h = 9(1 + 3^)(1-0)(1- 0/32 
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(l>4 = 9(1 + 0(1-»7)A/32 

4>5 = 9(1 +^)(1 - 377)(1 - 1;^)/32 
<i>^ = 9(1 + 0(1 + 3 » 7 )( 1 -,, 2)/32 

<i>7 = 9(1 + 0(1 +T?) a/32 
<!>% = 9(1+ 30(1 -e')(H- 0/32 
^9 = 9(1 - 30(1 -a(l + 0/32 
<^10 = 9(1-0(1+0A/32 
<kn = 9(1 -0(1 + 30(1 -t?')/32 
<I>12 = 9(1 -0(1 + 30(1 ->?')/32 

with 

A = (-10/9 + O + »7^)- 

Gauss Points and weights 

Gauss points 0 and their corresponding weights Wk used for numerical integration 
in the finite element method are given in Table B.l. The Gauss points and weights 
are given for different number of points N. 
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Table B.l; Gauss points and weights 


N 


Wk 

4 

-0.86113631159405 

-0.33998104358486 

0.33998104358486 

0.86113631159405 

0.34785484513745 

0.65214515486255 

0.65214515486255 

0.34785484513745 

5 

-0.90617984593866 

-0.53846931010568 

0.0 

0.53846931010568 

0.90617984593866 

0.23692688505619 

0.47862867049937 

0.56888888888889 

0.47862867049937 

0.23692688505619 

6 

-0.93246951420315 

-0.66120938646626 

-0.23861918608320 

0.23861918608320 

0.66120938646626 

0.93246951420315 

0.17132449237917 

0.36076157304814 

0.46791393457269 

0.46791393457269 

0.36076157304814 

0.17132449237917 
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